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1.0  INTRODUCTION 


The  influx  of  extra-terrestrial  material  is  responsible  for  the  dq)osition  of  some  sevraiteen 
thousand  tons  of  matter  in  Earth’s  atmosphere  each  year.  Among  the  important  consequorces 
of  this  process  is  the  formation  of  layers  of  free  m^s  at  around  100  km  altitude,  where  the 
atmosphere  becomes  dense  enough  for  substantial  ablation  of  the  particles  to  b^in.  The 
appearance  of  layers  of  neutral  atomic  metals,  «ich  as  sodium,  is  a  result  of  the  combination  of 
the  ablation  process  and  atmospheric  chemistry.  Ablation  produces  elem«ital  metal  atoms.  As 
the  metal  atoms  diffuse  downward,  they  are  increasingly  subject  to  three-body  reactions  as  the 
atmospheric  density  increases,  forming  the  various  metal  oxides.  From  diere,  a  chain  of 
reactions  results  evoitually  in  the  formation  of  conglomerates,  which  fall  to  Earth. 

Although  the  deposition  process  and  the  subsequoit  kinetics  of  the  metals  is  intimately  related, 
it  is  often  necessary  to  separate  the  modeling  effort  into  phases,  the  first  of  which  produces  a 
source  function  for  the  appearance  of  metal  atoms.  With  this,  one  can  next  address  the  diffusion 
problem  by  considering  the  transport  of  all  species  containing  a  particular  metal.  This  is  valid 
if  the  important  chemistry  takes  place  below  the  turbopause,  where  atmospheric  turbulence  leads 
to  a  mass  independent  s^e  height  for  minor  species.  Under  these  assumptions,  the  diffusion 
equation  can  be  solved  for  the  steady  state  concentration  of  the  total  number  of  atoms  containing 
the  metal.  Then,  a  kinetic  model  can  be  developed  and  solved  in  the  steady  state  to  give  the 
relative  amount  of  each  particular  molecular  metal  species.  Implicit  in  this  solution  is  the 
assumption  that  the  removal  of  all  metal  molecules  is  uniform,  which  is  not  completely  true  for 
the  complexes  that  are  the  eventual  result  of  the  chemistry.  The  assumption  is  quite  good  in  the 
metal  layer  itself,  however,  since  there  diffusion  is  r^id,  compared  to  the  non-diffusive  removal 
mechanisms. 

The  creation  of  a  model  can  be  divided  into  three  separate  problems:  the  generation  of  an 
altitude  dependant  source  function  giving  the  rate  of  dqx>sition  of  metal  atoms,  the  solution  of 
the  transport  equation  for  the  steady  state  total  number  density  of  metal  containing  species,  and 
the  creation  of  a  kinetic  model  to  give  the  relative  amounts  of  particular  species  at  a  particular 
altitude.  Different  modelers  have  in  the  past  approached  each  of  these  stq>s  in  rather  unique 
ways.  For  instance.  Plane  [1991],  in  his  model  of  sodium  deposition,  considers  eddy-diffusion 
only.  This  allows  the  entire  difliision  problem  to  be  based  on  only  an  assumed  flux  of  sodium 
containing  species  at  the  minimum  modeled  altitude  of  65  km.  This  flux  is  chosen  so  that  the 
model  results  in  a  reasonable  total  column  density  of  free  sodium,  matching  observations. 
Thomas,  et  al.  [1983]  began  with  a  source  profile  calculated  by  Hunten,  et  al.  [1980]  which  is 
scaled  uniformly  with  altitude  in  order  to  achieve  a  reasonable  free  sodium  peak  height.  Both 
the  treatments  neglect  molecular  diffusion,  which  limits  the  model  to  relatively  low  altitudes. 

In  the  present  work,  we  have  attempted  to  create  a  self-consistent  model,  beginning  with  the 
equations  describing  the  behavior  of  meteoric  particles  in  the  atmosphere.  By  following  the  time 
behavior  of  individual  particles,  we  are  able  to  derive  the  relative  amount  of  material  deposited 
by  them  at  various  altitudes.  This  is  then  coupled  with  the  estimated  size  distribution  functions 
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of  Hughes  [1992].  A  Monte  Cailo-type  calculation  is  tiien  performed  to  give  an  average 
dqxMiticm  as  a  function  of  altitude.  This  is  scaled  by  assuming  a  total  influx  of  material, 
estimated  by  several  woricers  to  be  around  507  grams  per  second  to  give  a  source  function.  The 
source  function,  obtained  in  total  grams  per  altitode  unit  per  sec<md,  is  then  i^lied  to  a  chosen 
metal  qxsdes  widi  known  average  weight  percents  of  die  chosen  metal. 

The  scniroe  function  is  then  combined  with  a  diffusion  model,  including  both  molecular  and  eddy 
diffusion,  to  give  steady  state  concentrations  of  conqiounds  of  die  chosen  metal.  A  kinetic 
model  is  dien  added  to  give  absolute  concentrations  of  various  metal  compounds.  We  test  this 
model  on  sodium  with  particular  attention  to  die  variation  in  sodium  ni^t-glow  with  ozone 
density.  The  chemistry  of  the  modd  is  introduced  in  bodi  a  simidistic  and  a  more  complicated 
sequence  of  reactions,  both  giving  results  in  good  agreement  with  obsovaticms. 


2.0  THE  DEPOSITION  FUNCTION 


The  first  stq>  in  the  modeling  is  the  calculation  of  the  rate  of  deposition  of  metal  atoms  at  a 
given  altitude.  The  rate  of  dqiosition  dqiends  on  the  atmospheric  density  and  the  vdocity  and 
size  of  the  particle.  To  a  lesser  extoit,  this  rate  dqiends  on  the  specific  compositicHi  of  the 
particle  as  well,  but  we  n^ect  this  complicaticm  for  simplidty.  The  particle  is  assumed  to  be 
spherical  in  shape  and  is  assumed  to  be  composed  of  a  material  like  Chmdrodite, 
[(Mg0H)2‘2Mg2Si04],  with  a  specific  gravity  of  3.2  g/cc.  The  velocity  of  the  particle  is 
subject  to  some  speculation.  If  the  particle  were  sitting  still  relative  to  the  when 
encountered,  the  Earth-fixed  velodty  would  be  almost  30  km/s.  But,  according  to  Hughes 
[1992],  most  particles  are  in  prograde  orbits  around  the  sun  (f.e. ,  they  revolve  around  the  sun 
in  the  same  sense  as  does  die  Earth).  This  makes  tiie  native  vdocity  substantially  less. 
Typical  geocentric  velodties  are  from  13.8  to  16.2  km/s,  again  according  to  Hughes  [1992]. 
We  will  investigate  the  velodty  dq)aidCTce  of  dqxmtion  and,  eventually,  find  that  the 
completed  modd  is  rather  insensitive  to  the  predse  vdodty  of  the  particle. 

The  evolution  of  a  particle  after  it  alters  Earth’s  atmo^ere  can  be  modeled  as  follows.  In  a 
time  df,  the  particle  oicounters  a  mass  of  atmo^ere  equal  to 

M„^A{m/pJ^p„Vdt  (1) 

where  is  the  atmospheric  density,  V  is  the  geocoitric  velodty,  m  is  the  mass  of  the  particle 
and  A  is  a  dimoisionless  shape  factor,  equal  to  about  1.2,  which  allows  the  expressicm  to  be 
writtoi  in  toms  of  mass  and  density  ratha  than  radius.  The  impact  of  the  particle  with  the  air 
serves  both  to  heat  the  air  and  the  particle,  leading  to  the  ablaticm  of  metal  atoms,  and  to  slow 
the  meteor  down.  The  change  in  velodty  with  time  contains  both  the  acceleration  of  the  particle 
due  to  gravity  and  the  atmospheric  drag. 
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where  tiie  first  term  is  tite  usual  expressicm  for  drag,  an  enqnrical  relationship,  and  the  second 
term  is  foe  acceleration  of  gravity.  Hughes  [1992]  gives  a  range  from  O.S  to  1.0  for  F  and  we 
choose  to  fix  the  value  at  0.7S  in  this  model.  The  accderation  due  to  gravity,  p,  is  4.0(20) 
cm^/s^. 

The  maximum  possible  energy  released  in  foe  air/particle  collision  is  equal  to  the  eneigy  that 
would  be  imparted  to  foe  air  molecules  if  they  were  raised  to  a  velocity  equal  to  the  particle 
velocity,  or, 


En»x  '  (3) 

which  follows  from  Eq(l).  It  is  assumed  that  some  fraction  of  this  enmgy,  A,  goes  toward 
heating  foe  material  in  foe  particle  instead  of  toward  heating  foe  air.  With  this,  one  can  write 
foe  mass  loss  due  to  ablaticm  as 


am 

■*"  2fpf 


where  |  is  the  heat  of  sublimation,  equal  to  around  5(10)  erg/g.  Hughes  [1992]  gives  a  range 
for  A  from  0.1  to  0.6.  We  fix  A  at  0.2. 

With  Eq^)  and  Eq(4),  along  with  the  time  dq)endrat  equation  for  altitude, 

^  =  V  (5) 

dt 

we  are  able  to  solve  for  foe  velocity,  altitude,  and  mass  of  a  particle  at  any  time,  givoi  its  initial 
velocity  and  mass.  The  only  elemoit  missing  is  foe  atmospheric  density  We  have  modeled 
this  af^  data  presented  by  Kelly  [1989]  for  sunspot  maximum.  Figure  1  shows  the  model 
density  in  g/cc. 

We  used  a  standard  4-fo  order  Runge-Kutta  integrator  to  solve  Eq(2),  Eq(4),  and  Eq(S) 
simultaneously.  Figure  2  foows  a  typical  plot  for  an  initial  mass  of  one  gram,  with  the  initial 
velocity  varied  between  10  and  20  l^s.  Above  about  90  km,  there  is  little  ablation  due  to  foe 
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Altitude  (km) 


low  atmospheric  d^isity.  The  bulk  of  the  ablation  takes  place  between  SO  and  80  km  altitude. 
Also,  below  about  SO  Im  altitude,  there  ^ypears  to  be  a  terminal  velocity.  The  mass  is  not,  in 
£u:t,  precisely  constant.  However,  the  vary  rapid  deceleration  which  takes  place  as  the 
atmo^here  b^ins  to  become  substantial  effectively  shuts  off  ablation  below  about  SO  km.  This 
also  greatly  reduces  the  particle  velocity,  which  would  explain  why  meteoric  particles  of  mass 
less  than  a  few  grams  are  only  luke-warm  when  they  reach  Earth. 

The  altitude  at  which  the  majority  of  mass  loss  takes  place  dq)ends  greatly  on  the  initial  mass, 
however,  the  ratio  of  terminal  mass  to  initial  mass  does  not  vary  greatly  with  terminal  mass. 
Figure  3  shows  the  mass  profiles  for  a  particle  originally  1(-S)  grams  in  mass.  The  greatest 
difference  between  this  and  the  1  gram  meteor  is  that  the  lighter  meteor  ablates  at  higho’ 
altitudes.  The  ablation  for  the  1(-S)  gram  particle  takes  place  between  about  70  and  110  km 
altitude.  The  mass  range  from  about  l(-8)  to  1  gram  is  the  most  important  r^on  for  this 
modeling,  since  this  range  contains  the  bulk  of  the  influx  material.  However,  it  is  interesting 
to  observe  this  behavior  in  extremes.  The  largest  meteorite  fall  in  the  British  Isles  was  the  48 
kg  meteorite  in  Limerick  in  1813  [Hughes,  1991].  Assuming  an  initial  velocity  of  10  km/s,  this 
meteor  would  have  ceased  to  ablate  at  around  10  km  altitude  and  would  have  slowed  to  about 
70  m/s  at  impact.  An  extremely  large  meteor,  on  the  other  hand,  would  not  decelerate 
appreciably  before  reaching  Earth’s  surface.  This  can  be  thought  of  as  an  extension  of  the 
behavior  shown  in  Figures  2  and  3,  in  which  the  region  of  ablation  falls  below  zero  altitude. 

The  variation  of  ablation  with  initial  velocity  is  also  evident  from  Figures  2  and  3.  There  are 
two  main  effects.  First,  the  fastest  particles  lose  nearly  all  their  mass  by  the  time  they  reach  the 
ground.  The  slowest,  however,  lose  only  about  75%  of  their  original  mass.  Also,  faster 
particles  ablate  at  slightly  higher  altitudes.  The  distribution  of  particle  velocities  is  actually  a 
very  poorly  known  quantity.  We  would  therefore  prefer  a  model  which  is  not  too  strongly 
d^)endent  on  the  initial  velocity.  We  will  be  able  to  factor  out  the  first  of  these  effects  in  the 
m^el.  The  precise  altitude  of  ablation,  however,  will  always  be  dependent  on  the  choice  for 
the  initial  velocity.  As  we  go  along,  we  will  often  present  results  for  initial  velocities  of  12,  14 
and  16  km/s,  presuming  that  these  rq)resent  minimum,  average  and  maximum  values.  If  we 
assume,  then,  that  the  calculation  can  be  carried  out  with  a  single  initial  velocity,  the  total 
deposition  can  be  modeled  from  the  distribution  of  initial  mass.  Fortunately,  this  distribution 
is  a  better  known  quantity  than  is  the  initial  velocity  distribution. 

The  calculation  of  the  total  deposition  at  a  specific  altitude  is  quite  straightforward,  if  we  know 
the  distribution  of  particle  mass  (frequoitly  called  0).  Hughes  [1975]  has  presented  a  detailed 
study  of  the  size  distribution  function  and  has  suggested  parameterizations  of  logjoO  that  are 
linear  in  log^^m.  We  have  carried  out  fits  to  his  data,  giving  the  following  functions  over  three 
mass  ranges. 


^ogio©  =  -  10.08  -  0.550  logiom 
logic®  =  ■  14-01  -  1-305  logjom 
logic®  =  -  15.00  -  1.500  logiom 
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10(-14)  <  m  <  10(-06) 
10(-06)  <  m  <  10(-02) 
10(-02)  <  m  <  10(2) 


Mass/Orighd  Mass 


Meteor  Deposition  for  M=1(— 5)  grams 


0.0  20.0  40.0  60.0  80.0  100.0  120.0  140.0 


Altitude  (km) 


Figure  3.  Same  as  Figure  2,  but  widi  an  initial  mass  of  grams. 
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6(m)  is  in  the  form  of  cumulative  number  flux  of  particles  greater  than  the  mass  m. 

Figure  4  shows  the  cumulative  flux  model,  the  units  of  which  are  particles-m‘^-s’^-2T-str'^  The 
actual  units  will  not  be  of  use  to  us,  since  we  will  normalize  to  the  total  influx  after  the 
calculation  of  rdative  dqx>sition.  Figure  S  shows  Ae  relative  number  of  particles  at  a  chosen 
mass,  multiplied  by  the  mass.  This  essentially  gives  the  cmitribution  of  a  particle  of  mass  m  to 
the  total  influx  of  mass.  Figure  5  shows  that  particles  betweoi  about  10(-7)  and  10(-2)  grams 
make  up  the  major  contributors  to  the  mass  influx  to  the  Earth. 

Now,  to  perform  the  calculation  of  deposition,  we  will  divide  altitude  into  a  series  of  bins,  each 
of  which  are  1  km  in  height.  Beginning  at  some  high  altitude,  we  will  select  representative 
particles  from  10(-14)  grams  to  10^)  grams  uniformly  represented  in  log^o  m.  We  will  follow 
each  of  these  re>resentative  particles  as  they  descend  through  the  atmosphere.  At  each  time 
step,  we  will  determine  which  altitude  bin  the  particle  is  in,  and  will  accumulate  within  that  bin 
the  mass  lost  in  the  present  timestep.  This  is  equal  to 

Am  =  AtW(m)  (6) 

at  . 

where  dm/dt  is  given  by  Eq(4)  and  W(m)  is  a  wdght  function  proportional  to  the  relative 
number  of  particles  of  a  particular  mass.  Numerically, 


W(m)  =  0(m)  -  0(m+fim) 

dm  is  a  computational  parametm’  represoiting  the  mass  st^  in  log^o  space.  We  choose  dm  to 
be  0. 1  for  the  most  detailed  calculations.  Along  with  the  altitude  bins,  we  also  accumulate  the 
total  mass  released  by  the  particle.  To  this,  we  add  the  mass  remaining  when  the  particle  hits 
the  ground  to  And  the  total  mass  dqposited  by  this  particle.  This  is  multiplied  by  the  weight 
function  and  summed  for  all  rqnesratative  particles  in  the  ensemble. 

Now,  according  to  several  audiors  [Hughes.  1975]  the  total  influx  to  the  Earth  of  inteiplanetaiy 
matter  is  around  507  grams/sec,  give  or  take  perhaps  150  grams/sec.  Therefore,  we  can 
normalize  the  total  influx  calculated  in  the  model  to  this  value  and,  if  we  normalize  the  altitude 
depoident  accumulations  by  the  same  factor,  the  result  will  be  the  deposition  rate  in  each 
in^vidual  altitude  bin.  The  validity  of  this  approach  is  pefliaps  best  shown  by  the  following 
reasoning.  The  accumulations  in  the  altitude  bins  give  the  amount  of  matter  deposited  in  a 
particular  bin  relative  to  the  total  matter  entering  the  atmosphere.  This  is  already  wdghted  for 
the  known  distribution  of  incoming  masses.  The  norm^ization  of  the  altitude  dependent 
depositions  and  total  mass  accumulated  in  the  calculation,  which  can  be  thougnt  of  as  taking 
place  ovCT  some  unspecifled  time  interval,  to  the  known  total  mass  influx  makes  the  altitude 
dqpendoit  results  equal  to  the  dqposition  rate. 
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From  there,  we  divide  by  the  volume  of  the  altitude  bin  to  obtain  the  grams  of  meteoric  material 
deposited  per  cc  at  individual  altitudes.  Next,  we  use  meteoric  composition  data  given  by  Plane 
[1991]  to  get  the  dqwsition  rate  in  atoms/cc.  The  wdght  percent  and  molecular  weights  of 
major  metallic  componoits  of  meteoric  material  are  given  in  Table  1. 


TABLE  1.  Meteoric  Metal  Composition 


Metal 

wt.  % 

MW 

Na 

0.6 

23.0 

Ca 

1.0 

40.1 

Ni 

1.5 

58.7 

A1 

1.7 

27.0 

Fe 

11.5 

55.8  1 

Mg 

12.5 

24.3  1 

This  concludes  the  calculation  of  the  deposition  rate  of  metals.  The  computation  itself  is  rather 
intensive,  requiring  about  four  hours  on  a  Sun  workstation.  The  time  required,  of  course, 
depends  on  the  choice  of  the  int^ration  time  stq)  and  the  mass  time  step  5m.  We  have 
investigated  the  dependence  of  the  calculations  on  both  these  parameters.  In  the  calculation,  the 
integration  timestep  is  chosra  such  that  thmre  will  be  200  Runge-Kutta  steps  for  each  Idlomet^ 
of  altitude.  This  depends,  of  course,  on  the  preset  velocity  of  the  particle.  The  mass  step  5m 
is  chosen  to  be  0.1.  This  gives  140  representative  particles  in  the  calculation  from  logjo^^  ~  “ 
14  to  log^om  =  2.  5m  was  chosen  small  oiough  so  that  the  diffmiences  in  ablation  height  were 
essentially  smoothed  out  over  the  1  km  altitude  bins,  leading  to  a  smooth  dqwsition  function 
with  altitude.  The  choice  of  5m  does  not  significantly  affect  the  deposition  rate  at  a  particular 
altitude,  as  long  as  the  deposition  rate  is  smooth  with  sdtitude,  since  the  dqx>sition  at  all  altitudes 
is  normalized  to  the  known  total  deposition  at  the  end  of  the  calculation.  A  final  parameter  in 
the  calculation  is  the  starting  altitude  of  the  particles.  We  have  found  that  starting  them  at  600 
km  is  high  enough  so  that  there  is  no  significant  change  in  the  results  that  would  have  been 
obtained  if  they  were  started  instead  at  500  km. 

The  deposition  function  will  be  called  q(h)  in  what  follows,  with  h  the  altitude.  Figure  6  shows 
q(h)  for  sodium  with  an  initial  particle  velocity  of  12  km/s.  Recall  that  we  do  not  attempt  to 
include  variations  in  initial  velocity  due  to  the  lack  of  knowledge  of  the  true  distribution. 
Instead,  we  set  the  initial  velocity  and  investigate  the  effects  of  changes  in  it  on  the  final  results 
of  the  model.  Figure  7  shows  a  blow-up  of  the  same  q(h).  The  shape  of  q(h)  is  really  quite 
similar  to  that  of  the  atmospheric  density  (see  Figure  1),  at  least  until  the  ablation  has  cea^, 
at  which  time  q(h)  rapidly  goes  to  zero.  The  peak  ablation  rate  for  Vo=  14  km/s  is  around  75 
km  altitude  with  a  rate  of  a  little  more  than  2(-2)  atoms/cm^-s.  Figure  8  shows  the  identical  plot 
for  Vo=12  km/s.  The  position  of  the  peak  is  little  changed  in  altitude,  but  the  deposition  rate 
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Altitude 


Meteor  Deposition  Model  with  V0=  14.0 
Compositbn  0.6  %  M.W.  23.0 


figure  6.  The  calculated  dqwsition  rate  iunction  for  sodium,  assuming  an  initial  particle 
velocity  of  14  km/s. 
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Altitude 


Meteor  Deposition  Model  with  VC)=  14.0 
Compositbn  0.6  %  M.W.  23.0 


Figure  7.  Same  as  Figure  6,  but  a  blow-up  of  the  region  of  substantial  Na  d^sition. 
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Meteor  Deposition  Model  with  V0=  12.0 
Compositbn  0.6  %  M.W.  23.0 


No  Deposition  Rote  (/cnn=i'*3-s) 


Figure  8.  Same  as  Figure  7,  but  with  an  initial  particle  velocity  of  12  km/s 


has  ^en  slightly,  by  about  S%.  We  will  not  dwell  on  the  compariscm  of  results  for  various 
Vg  values  at  this  point,  since  as  we  will  see,  q(h)  will  go  thrcHigh  two  further  int^radcms  before 
it  become  a  sodium  density.  It  will  be  more  expedient  to  wait  until  the  final  (tesired  quantity 
is  calculated  before  doing  further  comparisons.  Also,  we  will  not  introduce  the  results  for  rrther 
metals  here,  since  they  differ  from  so^um  only  by  a  scale  factor  which  can  be  calculated  easily 
from  the  data  in  Table  1. 

Before  moving  on  to  the  diffusion  problem,  we  will  perform  cme  more  calculation  widi  q(h). 
The  end  product  of  this  modd  is  a  steady  state  doisity  of  various  sodium  species,  NaX,  whoe 
X  rqnresaits  any  of  the  possible  sodium  spedes  present  in  the  atmosphere.  The  time  dq)aident 
diffiidon  problem 


I  .  (» 

where  w  is  the  transport  vdodty  and  n  the  concentration  of  NaX,  is  quite  difficult  to  solve  in 
goieral  and,  even  if  we  did,  we  would  still  want  to  carry  out  the  calculation  for  a  very  Itmg  time 
until  steady  state  doisities  were  obtained.  Thus,  we  seek  instead  a  solution  to  the  steady  state 
problem  from  the  outset.  In  order  to  remove  the  time  dq)endence  from  q(h),  we  notice  that  in 
the  steady  state,  the  flux  of  NaX  through  one  cm^  at  an  altitude  h  must  equal  the  rate  of  creation 
at  all  points  above  this  altitude.  If  this  were  not  so,  the  total  amount  of  NaX  above  any 
partici^  altitude  could  not  remain  constant.  In  fact,  since  Na  is  injected  at  all  altitudes,  it  must 
be  that  the  net  flux  of  NaX  is  downward  at  all  altitudes.  Otherwise,  a  steady  state  would  not 
be  possible.  Therefore,  we  will  concentrate  on  the  quantity  F(h),  the  downw^  flux  of  NaX. 

It  is  easy  to  calculate  F(h)  from  the  previous  result.  By  definition,  it  is  equal  to  die  creation  rate 
of  Na  at  all  points  above  the  altitude  h,  i.e.. 


(9) 


Figure  9  shows  F(h)  for  an  initial  particle  velocity  of  14  km/s.  As  expected,  the  downward  flux 
is  constant  below  about  60  km,  where  ablation  becomes  negligible.  Since  this  calculation 
requires  the  substitution  of  some  maximum  altitude  for  infinity  in  the  integral,  we  show  the 
results  of  variation  in  h^  in  Table  2.  Cleaiiy,  the  maximum  of  600  km  chosen  for  the  presort 
model  is  sufficient.  For  comparison.  Plane  [1991]  curtained  a  value  of  total  sodium  flux  of 
1.3(4)  cm'^-s'^  from  ablation.  His  result  was  obtained  by  adjusting  the  flux  at  65  km  solving 
upward  for  the  total  mixing  ratio  [see  Thomas,  1974]  making  the  total  column  density  of  sodium 
equal  presumed  experimoital  values.  We  obtain  a  total  downward  flux  (at  ground  level)  of 
alrout  4(4)  cm*^  s*S  which  is  about  a  factor  of  three  larger. 
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Altitude 


Meteor  Depoation  Mocte)  with  VC)=  14.0 
Compositbn  0.6  %  M.W.  25.0 


Downward  Rux  (/cm*=*2— s) 


Figure  9.  The  downward  flax  of  sodium  species  NaX,  at  various  altitudes  as  calculated 
from  die  deposition  model  with  =  14  km/s. 
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TABLE  2.  Variation  of  Downward  Flux  (/cm^-s) 
with  Int^ration  Cutoff 


Downward  Flux  at 

h,n.,=600  km 

hn^=5(X)  km 

h^=400  km 

1  km  altitude 

42769.29 

42769.29 

42768.98 

100  km  altitude 

943.0632 

943.9738 

943.6765 

Again,  as  with  q(h),  we  will  not  investigate  the  dependence  of  the  solution  on  the  choice  for 
initial  particle  velocity,  leaving  the  final  test  of  significance  until  we  obtain  chemical  species 
doisities.  We  are  now  ready  to  move  on  to  the  treatmoit  of  the  solution  of  the  steady  state 
diffusicHi  equation. 


3.0  DIFFUSION 


The  next  step  in  modeling  is  to  derive  from  the  source  function  a  steady  state  concentration  of 
NaX  at  various  altitudes.  By  solving  for  the  totid  concentration  of  Na  species,  we  avoid  at  this 
point  the  complications  of  the  kinetic  reactions.  This  a^^roach  to  the  problem  also  allows  us 
to  solve  for  only  one  minor  atmospheric  constituent,  NaX.  At  each  altitude,  the  source  of  NaX 
is  the  flux  from  altitudes  above  and  the  sink  doMmward  diffusion.  In  general,  the  bulk  motion 
of  the  NaX  depoids  on  pressure  gradients,  concentration  gradients,  thermal  gradients,  and 
differaitial  scale  heights.  When  a  species  is  a  minor  constituent  of  the  atmosphere,  as  is  NaX, 
the  diffusion  problem  can  be  simplitied  so  that  the  transport  velocity  w  can  be  written  in  terms 
of  the  concentration  gradient  of  NaX,  the  scale  height  of  the  bulk  of  the  atmosphere,  the  scale 
height  of  NaX,  and  the  temperature  gradient.  Following  Bcmks  and  Kockarts  [1973],  the 
expression  is 


w  =  -Di2 


1  ^  _1_ 
dh 


n 


\  dn^  ^  I  ^  I  dT 
Hi  dh  ^  Tdh 


(10) 


In  Eq(lO),  is  the  scale  height  of  NaX,  H  is  the  scale  height  of  the  atmosphere,  D}2  is  the 
molecular  dif^sion  coefficient,  K  is  the  Eddy  diffusion  coefficient,  nj  is  the  concentration  of 
NaX,  T  is  the  temperature,  and  h  is  the  altitude. 

The  mass  of  NaX  enters  the  equation  only  in  the  form  of  the  scale  height 
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(11) 


»1 


kr 


witii  g  tihe  aoceleraticm  of  gravity.  The  first  term  in  Eq(lO)  rqnesents  molecular  diffusirai, 
which  predominates  at  altitudes  of  above  about  1(X)  km.  The  second  term  rqnesents  Eddy 
diffusion,  predominant  in  lower  altitude  regions  ^>^iere  turbulence  tends  to  mix  all  spedes 
regardless  of  mass  .  We  can  get  away  widi  a  single  solution  to  this  equation,  using  only  the  scale 
height  for  atomic  sodium,  only  if  atomic  sodium  is  the  nugor  q)ecies  whenever  molecular 
diffurdm  predominates.  Althoti^  diis  turns  out  to  be  true  for  tee  sodium  model  presented  here, 
it  can  be  verified  only  in  retrospect. 

If  we  multiply  Eq(10)  by  tee  concentration  of  NaX,  nj,  we  obtain  a  very  ccmvenient  differential 
equation  n^, 


dn^ 

IR 


Dn 

(2>i2+X) 


K  r  1  ^  1  ^ 

(Di2+X)  [h  *  Tdhj  ^ 


m 

(I>12+X) 


i^teere  F  is  tee  downward  flux  function,  obtained  by  equating 

m)  »  "iW 


(12) 

(13) 


Eq(12)  is  instructive  since  tee  terms  in  [...]  are  really  just  scale  heights  which  contribute  to 
dn^/dh  in  relative  prc^xyrtion  to  tee  relative  magnitude  of  tee  Molecular  and  Eddy  diffusion 
coefficients  at  a  given  altitude.  In  tee  molecular  diffusion  regime,  tee  scale  height  is  teat  of  tee 
indqiendent  species,  while  in  tee  turbuloit  r^ime,  tee  scale  hdght  is  tee  same  as  that  of  tee 
general  atmo^teere. 

Another  notable  aspect  about  Eq(12)  is  tee  fact  teat,  should  q(h)  go  to  zero,  as  it  does  at  low 
altitudes,  tee  concentration  n^  still  increases  since  F(h)  goes  to  a  constant.  This  is  a  result  of 
tee  increasing  number  of  collisitnis  with  decreasing  altitude,  which  must  reduce  tee  transport 
vdodty  even  in  tee  absence  of  temperature  or  concentration  gradients.  Finally,  we  note  that, 
excqrt  for  tee  possible  thermal  ^adioit  term,  which  we  will  see  is  rdatively  small  by 
conpaiison,  tee  craicentration  n^  is  always  increasing  with  decreasing  altitude  according  to 
Eq(12). 

We  have  adopted  a  diffusion  coefficient  that  is  similar  to  that  of  argon  in  air,  as  rqx)rted  by 
Banks  and  Kockarts  [1973].  It  is  parameterized  by 
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(14) 


with  A  equal  to  5.5(16)  and  the  exponent  s  equal  to  0.841.  The  number  density  of  the 
atmosphere,  n,  is  calculated  from  the  model  mass  doisity  shown  in  Figure  1  and  the  altitude 
dq)aident  molecular  weight,  as  tabulated  by  Kelley  [1989].  Figure  10  shows  the  average 
molecular  weight  of  the  atmosphere.  Figure  11  shows  the  atmo^heric  temperature,  modeled 
from  the  same  source. 

Around  and  below  100  km  altitude,  turbulent  mixing  of  the  atmosphere  results  in  a  single  scale 
height  for  aU  species.  This  is  strictly  true  in  the  absence  of  chemical  reactions,  the  rates  for 
which  are  large  compared  to  the  eddy  diffusion  time.  For  this  model,  the  eddy  diffusion 
coefficient  was  modded  after  one  given  by  Thomas,  et  al.  [1983].  As  can  be  seen  in  Eq(12), 
the  diffusion  problem  dqrends  only  on  the  relative  magnitudes  of  the  molecular  and  eddy 
diffusion  coefficioits,  0^2  and  K,  and  on  the  sum  of  the  two,  compared  to  the  total  downward 
flux.  Figure  12  shows  the  relative  magnitudes,  and  Figure  13  shows  the  sum  of  the  two. 
Figure  11  is  limited  to  the  range  of  60  to  120  km,  since  molecular  diffusion  dominates 
completely  above  120  km,  and  eddy  diffusion  dominates  below  60  km.  This  r^on  is  also  the 
r^on  of  interest  for  the  chemical  model,  which  wUl  constitute  the  last  stq>  of  the  modeling 
effort.  It  is  important,  thoi,  that  both  molecular  and  eddy  diffusion  be  considmed,  since  the 
bulk  of  the  sodum  present  between  90  and  100  km  comes  from  above,  where  molecular 
diffusion  dominates.  Since  sodium  is,  in  general,  lighter  than  atmospheric  constituents,  it  is 
necessary  as  well  to  include  eddy  diffusion  to  oisure  a  net  downward  flux  of  sodium.  With 
molecular  diffusion  alone,  the  net  scale  height  of  the  sodium  would  become  greats  than  that  of 
the  bulk  atmosphere. 

The  atmospheric  scale  height  of  the  bulk  atmosph«e,  H,  is  givra  by 


H  =  0.871  ir/m  <15) 

with  H  in  km,  T  in  and  m  in  AMU.  ProfUes  for  T  and  m  have  already  been  presented.  The 
thermal  diffusion  term  is  calculated  by  numoical  differencing  of  the  temperature  profile.  The 
thermal  diffusion  fsictor  aj  is  takai  to  be  0.1,  about  that  for  Argon  in  air.  This  value  is  not 
critical  to  the  calculation  in  any  case. 

It  is  instructive  to  examine  the  relative  importance  of  the  terms  in  Eq(lO)  at  various  altitudes. 
To  do  so,  we  solve  the  equation  in  three  ways.  First,  we  use  only  the  source  term.  We  thai 
include  the  gravitational  terms  involving  the  scale  heights.  Finally,  we  solve  the  full  diffusion 
model  by  including  the  concentration  gradients.  Figure  14  shows  the  three  solutions  for  the 
initial  particle  velocity  14  km/s.  The  region  shown  is  from  50  to  2(X)  km  altitude,  whoe  the 
NaX  density  has  become  significant.  By  200  km,  we  see  that  the  density  of  solutions  2  and  3 
have  tripled  over  that  of  solution  1.  Solution  1,  with  the  source  t^m  alone,  represents 
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Atmospheric  Parameters  [after  Kelly,  1989] 


Figure  10.  The  average  molecular  weight  of  tbe  atmosphwe  [modeled  after  Kelley ^  1989]. 
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Atmospheric  Parameters  [after  Kelly,  1989] 


Tempercrture  (°C) 


Figure  11.  The  atmospheric  temperature,  in  "C  [after  Kelley,  1989]. 
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Molecular  and  Eddy  Diffusion  Terms 


Figure  12.  The  relative  magnitudes  of  die  eddy  and  molecular  diffiision  coefficients. 
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Total  Sodium  Density  Calculations 


10'"  v~'  rf  id  id  id  id  id  id 


NdX  Concentration  (/cc) 


Figure  14.  Steady  state  NaX  concentrations,  resulting  fromo  solution  of  Eq  (10)  widi  various 

terms  included  (see  text). 


2- 


something  of  a  ftee  molecukff  motion  where  the  NaX  is  not  influenced  by  its  own  concentradrm 
gradients.  The  solution  with  te]iq)erature  gradients  is  quite  similar  to  that  without  throughout 
the  altitude  range.  It  is  at  most  20%  higher  between  about  100  and  200  km,  where  the 
temperature  variation  is  greatest.  The  density  is  more  or  less  eqxxiential,  reaching  about 
3(2)/cc  at  100  km  and  5(3)  at  90  km,  the  expected  peak  of  the  sodium  layer. 

Anther  way  to  look  at  this  solution  is  to  compare  the  calculated  doisity  widi  the  ttMal 
atmosi^ieric  dmsity.  This  is  (kme  in  Figure  15.  We  see  that  the  mixing  ratio  is  shaped  more 
or  less  like  the  flux  functitm  in  Figure  9.  The  sdf-consistmcy  of  die  calculation  is  evident  in 
that  the  mixing  ratio  becomes  constant  below  the  region  of  dep^tion.  This  means  that  the  NaX 
distribution  has  adsqited  die  scale  height  of  the  general  atmosphere,  and  the  source  term  in 
Eq(12)  has  become  insignificant  by  comparison.  The  ^‘roughness”  in  the  low  altitude  r^cm 
of  Figure  15  is  a  result  of  the  very  h^e  dentities  there  and  the  discreetness  of  the  total 
atmosphere  model.  We  do  not  know  if  1(-10)  is  a  reasonable  mixing  ratio  for  sodium  species 
at  ground  level,  but  at  least  it  does  not  seem  Mcessively  large.  We  also  note  that  there  are 
heterogoieous  removal  mechanisms  in  the  troposphm  which  no  doubt  can  remove  sodium 
species  far  tiister  than  diffusion  below  50  km  or  so. 

A  third  way  to  examine  this  solution  is  to  compare  the  scale  height  of  the  sodium  to  that  of  the 
total  atmosphere.  Figure  16  shows  these.  At  first,  with  the  calculation  beginning  at  600  km, 
the  scale  height  is  low.  It  increases  until  about  350  km,  at  which  point  diffusicm  b^ins  to 
dominate  the  creation  from  ablation.  From  350  km  to  about  80  km,  the  scale  height  of  sodium 
decreases  along  with  the  total  atmosphere  scale  height.  Below  80  km,  it  adopts  £q>proximately 
the  same  scale  hdght  of  the  atmosphere,  due  to  eddy  diffusion.  For  this  rough  comparison,  the 
scale  height  was  computed  numerically  from  successive  points  by  the  formula 


=  A/t/Alogn(A/i) 

Since  this  calculation  dq)aids  on  the  stq)  size  and  the  starting  altitude,  we  should  check  it  with 
various  values  of  step  size  Ah  and  start  altitude  Iiq.  Table  3  shows  a  comparison  of  the  result 
that  will  turn  out  to  be  critical  to  the  kinetic  model,  the  NaX  doisity  at  90  km,  for  various  step 
sizes  and  initial  altitude.  In  Table  2,  NaX  doisities  are  in  fen?,  and  stq>  size  and  start  altitude 
in  km.  Clearly,  considering  the  precision  of  this  model,  the  integration  parameters  are  adequate. 


TABLE  3.  NaX  Density  with  Various  Inte 

gration  Parameters 

Step  Size  — > 

1  km 

0.2  km 

ho  =  6(X)  km 

2.655(3) 

2.667(3) 

ho  =  500  km 

2.652(3) 

2.660(3) 

ho  =  400  km 

2.642(3) 

2.646(3) 
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Mitude  (km) 


NdX  Diffuabn  Results 


Mixhg  Ratio 


Figure  15.  Calculated  mixing  ratio  of  total  NaX  density  from  solution  of  die  diff-.!^.^ 
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Altitude  (km) 


Having  now  arrived  at  a  steady  state  density  profile  of  total  sodium  q)ecies  as  a  function  of 
altitude,  we  are  ready  to  move  on  to  incorporate  the  kinetic  reactions  of  sodium  with 
atmoq^eric  cmistituents. 

4.0  THE  KlNEnC  MODEL 


The  first  kinetic  modd  examined  is  a  simplified  out,  similar  to  that  used  by  Smder  [1986]  in 
his  study  of  sodium  nightglow.  This  kinetic  modd  contains  mily  dght  reacticms  between  four 
sodium  species  and  the  ambient  atoms,  but  it  contains  all  the  essential  physics  necessary  to 
obtain  reasonable  estimates  of  the  sodium  layer  and  die  nightglow  intensity.  It  is  also  useful  in 
understanding  the  general  features  of  the  sodium  layer  because  of  its  relative  simplidty.  When 
we  move  on  to  a  more  conqilicated  modd,  we  will  see  that  the  results  from  this  simplm*  one  are 
particularly  useful  in  adjusting  the  rate  omstants,  since  the  simpler  model  is  far  less  sensitive 
to  small  rate  changes  than  is  the  more  complex  one. 

The  origin  of  the  S894A  nightglow  was  postulated  by  Chapman  [1939]  as  due  to  the  reaction 
between  NaO  and  O  to  form  atomic  sodium,  a  fraction  of  which  is  left  in  the  excited 
electronic  state. 

NaO  +  O  Na(2p,2s)  +  O2 

The  exdted  state  would  radiate  almost  immediately  to  the  ground  state.  Cluq>man  further 
proposed  that  NaO  would  be  generated  from  atomic  Na  via  the  reaction 

Na  +  03**  NaO  +  O2  , 

another  mechanism  for  formation  of  NaO  is  the  three  body  reaction 

Na  +  O  +  M  -♦  NaO  +  M  . 


But  this  reaction  is  of  secondary  importance  in  the  process,  as  discovered  by  Bates  and  Ojha 
[1980].  Hussin  and  Plane  [19S2]  later  suggested  a  diird  mechanism  involving  the  reaction 

Na  +  O2  +  N2  Na02  +  N2 


followed  by 


Na02  +  O  -*  NaO  +  O2 


for  NaO  formation. 
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The  five  reactiims  above  would  constitute  a  complete  model  for  the  steady  state  behavior  of 
sodium.  However,  if  we  went  ahead  and  solved  this  modd,  we  would  find  that  the  profile  of 
free  sodium  does  not  at  all  resemble  a  layer,  as  is  found  experimentally.  What  is  laddng  so  far 
is  a  reaction  to  turn  sodium  into  a  relatively  long  lived  q)ecies  at  Iowa-  altitudes.  The  complex 
chain  of  reactions  most  likely  involves  more  and  mme  stable  qiecies  as  the  reactions  progress, 
but  one  stable  q)edes  that  can  be  created  directly  is  NaOH.  This  is  created  from  NaO  by  die 
reaction 


NaO  +  HjO  -  NaOH  +  OH  , 
and  from  Na02  by  the  reaction 

NaOj  +  OH  NaOH  +  Oj  . 

These  reactions  have  reasonably  well  known  rate  constants.  They  affect  the  Iowa*  altitude 
densities  of  free  sodium  because  both  OH  and  H2O  increase  with  decreasing  altitude.  Finally, 
in  order  to  solve  the  model,  we  need  to  add  a  ‘closure’  reaction.  If  we  did  not,  then  the  steady 
state  solution  would  result  in  100%  NaOH,  since  NaOH  is  never  returned  to  any  of  the  other 
three  constitu^its.  A  reasonable  reaction  for  this  purpose  is 

NaOH  +  H  Na  +  HjO  . 

Although  the  rate  of  this  reaction  is  not  well  determined  (and  in  fact,  other  reactions  producing 
more  complex  sodium  compounds  may  dominate  it),  this  does  not  matter,  since  we  will  show 
that  the  characteristics  of  the  sodium  layer  and  the  nightglow  are  not  strongly  dq)aident  on  this 
reaction  rate.  It  is  necessary  only  for  closure.  The  reactions,  and  rate  constants  as  suggested 
by  Swider  [1986]  are  given  in  Ts^le  4. 


TABLE  4.  Kinetic  Model  A 


B 

Reaction 

Rate  1 

1 

Na  +  O3  -*  NaO  +  O2 

7(-10) 

2 

NaO  +  0  -  Na  +  O2 

K-IO) 

3 

Na  +  O2  +  N2  Na02  +  N2 

2(-30) 

4 

NaO  +  H2O  -  NaOH  +  OH 

2(-10) 

5 

NaO  +  O3  -*  Na02  +  O2 

2(-10) 

6 

Na02  +  0  NaO  +  O2 

1(-13) 

1  ^ 

Na02  +  OH  NaOH  +  O2 

K-IO) 

1  S 

NaOH  +  H  Na  +  H2O 

K-ll) 
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In  order  to  proceed,  we  must  make  up  altitude  dq)aident  models  of  the  major  and  minor 
constitumts  of  the  ambient  atmosphere.  To  simplify  things,  we  will  consider  only  the  region 
betwera  60  km  and  120  km  altitude.  Since  the  steady  state  problem  has  now  been  reduced  to 
a  total  density  for  NaX  at  fixed  altitude,  the  integrations  over  altitude  are  complete,  and  the 
problem  is  now  solved  independently  at  each  altitude.  In  what  follows,  we  will  need 
conc«itrations  of  N2,  O2,  H,  OH,  H2O,  O,  O3,  CO2,  02'*',  NO'*',  and  e'.  Although  we  do  not 
need  them  all  for  the  present  model,  we  will  present  them  all  here  for  simplicity. 

First,  we  note  that  we  are  most  interested  in  calculation  of  the  Na(D)  nightglow.  Thus,  v/e  will 
restrict  the  modeling  to  the  night  time.  For  the  major  atmospheric  constituents  N2  and  O2,  we 
assumed  perfect  mixing  below  83  km  altitudes,  with  the  U.S.  Standard  Atmosphere  mixing  ratios 
of  78.06%  and  20.94% ,  respectively.  Above  83  km,  the  densities  tabulated  in  the  U.S.  Standard 
Atmosphere  were  used  directly.  These  were  interpolated  linearly  in  log^onfh).  Figure  17  shows 
the  result.  Minor  constituents  H,  O,  and  OH  were  modeled  after  polynomial  fits  to  discrete 
points  taken  from  the  results  of  modeling  of  Alkn,  et  al.  [1984].  Figure  18  shows  these 
concratrations.  For  H2O,  the  mixing  ratio  reported  by  Allen,  et  al.  was  fit  to  a  polynomial  and 
then  multiplied  by  the  total  atmospheric  density,  which  was  in  turn  obtained  from  the  model 
atmo^here  in  Figure  1  and  the  average  molecular  weight  in  Figure  10.  The  CO2  profile  was 
produced  by  another  polynomial  fit  to  the  curve  presented  by  Keneshea,  et  al.  [1979].  Figure 
19  shows  these  d^sities. 

The  ozone  density  was  allowed  to  vary  somewhat  in  these  models,  since  part  of  the  purpose  of 
the  kinetic  modeling  was  to  ascertain  the  effect  of  varying  ozone  density  on  the  sodium  night- 
glow.  We  adapted  a  general  profile  reported  by  Allen,  et  al.  [1984]  as  determined  from  rocket 
measuremeits,  but  allowed  the  density  to  vary  by  a  factor  of  10  over  that  reported.  Figure  20 
shows  the  minimum  and  maximum  ozone  conceitrations  considered  in  this  model.  For 
reference,  our  minimum  ozone  density  at  90  km  is  about  a  factor  of  two  less  than  that  of  Sze, 
et  al.  [1982].  Also,  McPeters  [1980]  has  noted  an  approximately  10-fold  increase  in  O3  at  1 
mbar  (altitude  —  50  km)  from  Nimbus-4  BUV  measuremoits,  so  this  range  of  variability  is  not 
eqiecially  extreme.  Finally,  we  will  need  to  model  charge  exchange  and  molecular  dissociative 
recombination,  for  which  we  take  NC*",  02''',  and  e*  densities  from  the  International  Reference 
Ionosphere  (IRI)  at  conditions  for  January  at  midnight  at  the  equator.  These  are  shown  in 
Figure  21. 

To  solve  this  model,  and  all  subsequent  ones,  a  file  is  prepared  containing  all  reactions  and  rate 
constants  to  be  considered.  We  have  written  a  software  program  called  KINO  which  reads  this 
file  and,  from  it,  prq)ares  a  subroutine  through  which  the  equations  are  solved.  This  prevents 
errors  in  coding,  which  can  become  quite  complicated  as  the  number  of  reactions  and 
constituents  increase.  The  system  of  equations  is  Aen  solved  in  the  steady  state.  The  details 
of  the  KINO  routine  and  of  the  steady  state  solution  are  givra  in  the  appendix. 

The  steady  state  solution  to  Model  A  is  shown  in  Figure  22.  The  kinetics  has  produced  a  well 
defined  free  sodium  layer,  centered  around  about  88  km  altitude.  Below  this,  there  is  a 
substantial  concmtration  of  both  Na02  and  NaOH.  At  even  lower  altitudes,  the  Na02  has  given 
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Night-time  Modd  Atmosphere 


Figure  17.  N2  and  O2  3S  a  function  of  altitude  for  tiie  night-time  model  atmosph^e. 
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Figure  18.  O,  H,  and  OH  doisity  as  a  function  of  altitude  for  the  night-time  model  atmosphere. 
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Figure  20. 
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Minimum  and  maximum  modded  O3  density  as  a  fimetion  of  altitude  for  foe 
ni^t-time  model  atmosphere. 
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Figure  21.  Model  ionosphere  OKI  at  midni|^t,  equator,  January)  as  a  function  of  altitude. 
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way  in  favor  of  NaOH  completely.  The  reason  that  there  is  a  layer  of  free  sodium  is  that  atomic 
oxygen  rayndly  converts  NaO  back  to  free  sodium  and  atomic  hydrogoi  NaOH  to  Na  at  these 
altitudes.  As  the  concentration  of  O  and  H  dr(q)s  off  below  85  km  or  so,  these  reactions  no 
longer  dominate.  NaOH  becomes  the  major  product  at  low  altitudes  because  it  is  destroyed  only 
by  the  reaction  with  H,  which  trails  off  even  more  rapidly  than  O  or  OH,  according  to  Figure 
18.  NaOH  is  the  major  low  altitude  species  in  this  model  rnily  because  no  further  reactions  have 
been  specified  which  could  convert  it  into  more  complex  ^)ecies,  NaHCOs  for  example.  Evoi 
with  its  simplicity,  though,  the  model  has  rq)roduced  the  sodium  layer  well,  in  agreemoit  with 
experimoital  evidence.  Tilgner  and  von  Zoito  [1988],  for  example,  have  calculated  the  mean 
«xlium  layer  at  64**  north  in  winter  from  lidar  measurements.  Their  result  shows  a  peak  height 
of  around  88  km  and  a  peak  magnitude  of  from  3(3)  to  4(3)  particles/cm^.  Model  A,  as 
evaluated  here  with  an  initial  particle  velocity  of  14  km/s,  gives  a  peak  hdght  at  about  87  km 
with  magnitude  of  3.3(3)  particles/cm^. 

Here,  we  investigate  the  behavior  of  the  model  as  a  function  of  particle  velocity.  As  we  recall, 
the  initial  velocity  of  the  particle  was  held  fixed  in  the  Monte  Carlo  calculation  of  deposition. 
This  was  mainly  because  no  data  was  available  on  the  precise  distribution  of  these  velocities. 
Hughes  [1991]  states  that  a  reasonable  range  for  the  vast  majority  of  particles  is  13.8  to  16.2 
km/s.  In  Table  S  we  show  the  various  values  of  the  peak  intensity  of  the  sodium  layer  for  initial 
velocities  between  12  and  16  km/s.  These  are  shown  both  for  ozone  minimum  and  ozone 
maximum  values.  In  all  cases,  the  peak  height  is  dther  87  or  88  km  altitude.  As  can  be  seoi, 
the  dependoice  on  initial  velocity  is  small,  less  than  10%  in  all  cases.  Intoestingly  too,  the 
maximum  in  the  sodium  layer  intensity  comes  at  around  IS  km/s.  The  variation  exhibits  a 
maximum  because  the  ablation  height  falls  as  velocity  increases,  however,  as  the  height  falls, 
more  Na  becomes  Na02  and  NaOH,  leading  to  an  overall  decrease.  On  the  other  side,  ablation 
at  higher  altitudes  increases  peak  density. 


TABLE  5.  Variation  of  Na  Layer  Maximum  (/cc) 
with  Initial  Particle  Velocity 


Velocity 

Ozone 

Ozone 

(km/s) 

Minimum 

Maximum 

12 

3.037(3) 

2.104(3) 

14 

3.240(3) 

2.244(3) 

16 

3.226(3) 

2.233(3) 

Considering  now  the  problem  of  sodium  night-glow,  we  note  that  the  number  of  photons  per 
cm^-s  should  be  proportional  to  the  rate  of  reaction  2  in  this  model.  This  proportionality 
constant  will  be  the  branching  ratio  between  the  and  states.  It  is  well  known  that  about 
one-third  of  the  Na  from  reaction  2  ends  up  in  the  excited  ^P  state,  so  that  the  volume  emission 
rate  is  given  by 
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Figure  23  shows  the  volume  emission  rate  as  a  function  of  altitude  for  this  model  at  ozone 
minimum.  The  night-glow  layer  peaks  at  about  35  {rfiotons/cm^-s  at  a  hdght  about  equal  to  die 
height  of  the  sodium  layer,  88  km. 

To  arrive  at  a  total  intensity  of  Na(D)  onission,  we  simply  int^rate  up  the  column. 

OD 

7=1  rih)dh 


which  we  do  by  Simpscm’s  rule.  Since  one  Rayldgh  is  1(6)  photons/cm^-s  and  since  we  use  1- 
km  altitude  divisions,  the  normalization  factor  is  one-tenth.  For  the  case  shown  in  Figure  23, 
the  total  intensity  is  24.7  R,  for  ozone  minimum.  In  the  case  of  ozone  maximum,  the  intoisity 
comes  out  to  be  130.2  R,  an  increase  of  about  500%.  The  increase  in  night-glow  wiA 
increasing  ozone  ccmcentration  is  not  direcdy  prqxirtional  because,  as  we  can  see  from  Table 
5,  the  increase  in  ozone  causes  a  decrease  in  steady  state  Na  density,  which  in  turn  drives  the 
emission  process.  Also,  Reaction  3  of  Model  A,  followed  by  Reaction  6  is  also  a  source  of 
NaO,  although  a  relatively  minor  one.  We  can  examine  this  by  running  the  model  without  any 
ozone  at  all.  The  major  effects  are  the  lowering  of  the  night-glow  peak  to  about  85  km  and  the 
lowering  of  the  intensity  to  about  4.1  R.  Thus,  at  ozone  minimum,  the  three-body  reaction 
process  contributes  about  20%  to  the  total  night-glow. 

In  the  process  of  creating  a  more  comprdiaisive  model  of  sodium  chemistry,  we  will  add 
reactions  to  the  chain  a  few  at  a  time.  We  now  introduce  Model  B,  which  will  add  three 
species.  First,  we  will  attempt  to  rqnoduce  the  model  of  Plane  [1991],  who  includes  a  total  of 
sevoiteen  reactions.  These  are  shown  in  Table  6.  Aside  from  a  few  minor  changes  in  rale 
constants,  there  are  some  differences  between  the  reactions  in  Model  A.  The  reaction  of  NaO 
with  ozone 


NaO  +  O3  -*  Na  -I-  2O2 

has  been  added  as  a  way  to  recycle  atomic  sodium.  Also,  three-body  reactions  betwera  NaO 
and  O2  or  CO^  have  bem  added,  and  a  three-body  reaction  between  NaOH  and  CO2,  as  well. 
Also  of  note,  the  recycling  of  Na02  with  OH  has  been  rqilaced  by  recycling  with  H  in  Model 
B.  The  three-body  reactions  of  NaO  now  form  NaO^  and  NaC03.  Both  of  these  are  convoted 
back  to  Na02  with  atomic  oxygen.  Finally,  NaCQj  is  allowed  to  react  with  hydrogen,  also  in 
a  three-body  reaction,  to  form  NaHC03.  NaHC03  is  recycled  to  atomic  sodium  with  atomic 
hydrogen  alone. 

NaHC03  +  H  -  Na  +  H2CO3 
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Figure  23. 


Knetic  Model  A  —  O3  Minirmum 


The  rate  of  Na(^P)  emission  as  a  fimction  of  altitude,  calculated  from  die  solution 
to  Kinetic  Modd  A  at  O3  minimum- 


We  will  see  that  this  reactim  is  especially  important  in  the  overall  behavior  of  this  model. 


Figure  24  shows  the  sodium  q)ecies  resulting  from  the  solution  of  Model  B.  The  first  striking 
difference  betwemi  Model  B  and  Model  A  is  that  the  peak  sodium  density  in  Model  B  has 
decreased  to  8.2(2)/cm?.  The  height  of  the  layer  has  increased  as  wdl  to  93  km.  At  90  km, 
the  major  species  is  NaHCOs,  comprising  more  than  80%  of  the  total  NaX.  The  intensity  of 
the  Na(D)  emissicm  has  falloi  dramatically  to  only  2.2  R. 


TABLE  6.  Kinetic  Model  B 


# 

Reaction 

Rate 

1 

Na  +  O3  -•  NaO  +  O2 

6.0(-10) 

2 

NaO  +  0  Na(2p,2s)  +  0^ 

2.2(-10) 

1  ^ 

NaO  +  03-*  Na02  +  O2 

1.5(-10) 

1  ^ 

NaO  +  O3  -*  Na  +  2O2 

2.0(-ll) 

1  S 

NaO  +  H2O  -*  NaOH  +  OH 

1.8(-10) 

6 

NaO  +  O2  +  N2  Na03  +  N2 

5.3(-30) 

7 

NaO  +  CO2  +  N2  NaC03  +  N2 

1.3(-27) 

8 

Na  +  O2  +  N2  -  Na02  +  N2 

4.7(-30) 

9 

Na02  +  0  NaO  +  O2 

2.0(-14) 

10 

Na02  +  H  -*  NaOH  +  0 

1.9(-14) 

11 

NaOH  +  H  -  Na  +  H2O 

4.0(-13) 

12 

NaOH  +  CO2  +  N2  NaHC03 

1.9(-28) 

13 

Na03  +  0  Na  +  2O2 

2.0(-14) 

14 

Na03  +  0  -*  Na02  +  O2 

2.0(-14) 

15 

NaC03  +  0  Na02  +  CO2 

1.0(-13) 

16 

NaCOa  +  H  +  N2  -  NaHC03  +  N2 

1.0(-30) 

17 

NaHC03  +  H  -  Na  +  H2CO3 

1.0(-14) 

Comparing  these  results  to  those  of  Plane,  we  find  that  the  peak  height  agrees  quite  well. 
However,  he  obtains  a  peak  density  of  about  3(3)/cm^  and,  more  importantly,  the  sodium  at  the 
peak  is  almost  entirdy  in  the  atomic  form.  His  peak  NaO  density  is  approximately  10/cm^, 
leading  to  an  integrated  intensity  of  106  R,  as  reported  there.  Finally,  the  model  evaluated  with 
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OUT  atmosi^iere  shows  that  the  major  constituoit  below  the  sodium  layer  is  NaHC03.  Howevo', 
the  figure  in  Plane  seems  to  indicate  that  Na02  is  tl»  major  species.  We  believe  that  there  must 
be  some  dramatic  difference  in  our  atmospheric  modd  from  that  used  by  Plane.  This  is 
somewhat  difficult  to  verify,  since  the  atmosphere  used  by  Plane  is  not  explicitly  specified. 

It  is  instructive  to  examine  the  effects  of  variaticms  of  cmlain  key  rate  cmistants  in  the  model, 
in  order  to  ascertain  the  reasons  for  the  possible  discrq>ancies.  First,  we  note  that  reaction  9, 
between  Na02  and  O  is  a  factor  of  five  smaller  in  Model  B  than  in  Model  A.  We  try  increasing 
it  by  a  factor  of  five.  The  effect  is  shown  in  Table  7.  This  reaction  increases  the  peak  density 
and  radiance  only  slightly.  Next,  we  examine  the  ‘‘closure”  reactions,  for  which  the  rates  are 
very  speculative  for  lack  of  measurements.  Increasing  Reaction  IS  and  Reaction  16  rates  by 
factors  of  100  has  almost  no  effect.  However,  increasing  the  rate  at  which  NaHCO^  is 
converted  back  to  Na  increases  the  peak  daisity  to  almost  what  it  was  in  Model  A.  Importantly, 
also,  we  see  that  after  the  rate  is  increased  1,000  times,  further  increases,  evoi  by  a  factor  of 
10,  do  not  significantly  change  the  results.  This  is  as  it  should  be,  indicating  that  atomic  Na  is 
the  major  species  in  the  lower  thermosphere. 


TABLE  7.  Effects  of  Variation  of  Rates  in  Model  B 


1  Variation 

Na  Peak  (/cc) 

Peak  Altitude 
(km) 

Night-glow 

(R) 

Original  Model 

818 

93 

2.4 

Reaction  9  increased  by  a 
factor  of  5 

959 

93 

2.8 

Reaction  IS  increased  by  a 
factor  of  100 

959 

93 

2.8 

Reaction  16  increased  by  a 
factor  of  100 

959 

93 

2.8 

Reaction  17  increased  by  a 
factor  of  100 

2313 

89 

16.3 

Reaction  17  increased  by  a 
factor  of  1,000 

2588 

88 

22.2 

Reaction  17  increased  by  a 
factor  of  10,000 

2630 

88 

23.5 

We  cannot  determine  precisely  why  the  Plane  kinetic  model  does  not  produce  results  in  good 
agreement  with  the  simpler  Model  A.  However,  the  adjustment  of  the  closure  rates  does  bring 
the  results  close  in  line  with  both  the  simpler  model  and  with  expected  results  based  on 
observation. 
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TABLE  8.  Kinetic  Model  C 


# 

Reaction 

Rate 
Model  B 

Rate 

Model  C 

1 

Na  +  O3  NaO  +  O2 

6.0(-10) 

7.0(-10) 

2 

NaO  -t-  0  -♦  Na  +  O2 

2.2(-10) 

1.0(-10) 

3 

NaO  +  O3  Na02  +  O2 

1.5(-10) 

2.0(-10) 

4 

NaO  -1-  O3  -*  Na  2O2 

2.0(-ll) 

2.0(-ll) 

5 

NaO  +  HjO  -*  NaOH  +  OH 

1.8(-10) 

2.0(-10) 

6 

NaO  -f-  O2  +  N2  Na03  +  N2 

5.3(-30) 

5.3(-30) 

7 

NaO  +  CO2  +  N2  NaC03  +  N2 

1.3(-27) 

1.3(-27) 

8 

Na  +  O2  +  N2  Na02  +  N2 

4.7(-30) 

2.0(-30) 

9 

Na02  +  0  NaO  +  O2 

2.0(-14) 

1.0(-13) 

10 

Na02  +  H  -*  NaOH  -1-  O 

1.9(-14) 

1.9(-14) 

11 

NaOH  +  H  -*  Na  +  H2O 

4.0(-13) 

1.0(-11) 

12 

NaOH  +  CO2  +  N2  -*  NaHCOj 

1.9(-28) 

1.9(-28) 

13 

Na02  +  OH  -*  NaOH  +  O2 

— 

1.0(-10) 

14 

Na03  +  0  -»  Na  +  2O2 

2.0(-14) 

2.0(-14) 

15 

Na03  -!■  0  -*  Na02  +  O2 

2.0(-14) 

2.0(-14) 

16 

NaC03  -HO-* Na02  +  CO2 

1.0(-13) 

1.0(-11) 

17 

NaC03  -H  H  -H  N2  -*  NaHC03  +  N2 

1.0(-30) 

1.0(-28) 

18 

NaHC03  -H  H  -*  Na  -H  H2CO3 

1.0(-14) 

1.0(-11) 

We  make  a  few  more  minor  modifications  to  this  model  to  arrive  at  Model  C,  which  is  shown 
in  Table  8.  The  modifications  include  rate  constant  modifications  to  conform  better  to  Model 
A,  and  the  re-introduction  of  the  reaction 


Na02  -I-  OH  NaOH  -I-  O2  , 

which  was  in  Model  A  but  not  in  Model  B.  For  comparison.  Table  vm  also  shows  the  rate 
constants  from  Model  B  alongside.  The  results  for  the  constituents  from  the  solution  of  Kinetic 
Model  C  are  shown  in  Figure  25.  Notably,  NaHC03  remains  the  predominant  species  below 
the  sodium  layer.  NaO,  which  is  responsible  for  the  sodium  nightglow,  has  come  back  up  to 
a  level  of  about  10/cc  at  an  altitude  of  85  km. 
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In  the  next  step,  we  examine  further  what  might  ha|^)ai  to  the  primary  low  altitude  product 
NaHC03.  NaO^  reacts  with  O,  H,  and  OH  at  high  altitudes.  However,  at  lower  altitudes  there 
is  very  little  of  these  three  species  present.  There,  the  primary  reaction  is  probably  the  three- 
body  reaction 


Na02  +  CO2  +  N2  NaC03  -I-  O  +  N2  , 
and  NaC03  can  react  further  through  another  three-body  reaction 

NaC03  +  Na  +  N2  -*  Na2C03  +  N2  . 

Eventually,  the  Na2C03  would  form  complex  crystals  and  ^  from  the  atmosjrfiere.  However, 
the  modeling  of  this  is  somewhat  complex  and  we  will  terminate  the  chain  of  reactions  here. 
In  order  to  solve  this  model,  however,  we  need  a  closure  reaction,  which  we  choose  to  be 

Na2C03  -t-  O  NaC03  +  NaO 

Choosing  rate  constants  of  1.0(-30),  1.0('‘28),  and  1.0(-16)  for  these  three  new  reactions 
respectively,  gives  the  results  for  our  Kinetic  Model  D,  shown  in  Figure  26.  Although  the 
results  above  about  80  km  have  not  changed  at  all  as  a  result  of  the  addition  of  the  three-body 
reactions,  we  see  that  the  low  altitude  model  is  quite  different.  The  NaHC03  has  formed  a  layer 
with  a  maximum  of  about  l(5)/cc  at  around  70  km.  Below  this,  the  major  component  is  NaC03 
for  a  short  time,  giving  way  at  lower  altitudes  to  Na2C03.  Increasing  the  rate  of  the  second  of 
the  new  reactions  in  relation  to  the  first  would  no  doubt  reduce  the  NaCO^.  However,  without 
either  accurate  measurements  of  the  rates  or  experimental  observations  of  these  species,  this 
would  be  a  rather  speculative  exercise. 


Finally,  we  finish  out  the  model  by  adding  the  ionic  reactions.  These  include  charge  exchange 
of  atomic  Na  with  both  02'‘'  and  NO'*', 


and 


Na  +  02^  Na-*-  +  O2 

Na  +  NO-*-  -  Na-*-  +  NO  . 


The  sodium  ion  is  then  assumed  to  react  with  two  N2  atoms  to  form  the  relatively  stable  complex 
NaN2''' 


Na-*-  -1-  N2  +  N2  -  NaN2'''  +  N2  . 

Finally,  molecular  dissociative  recombination  (MDR)  can  take  place  to  close  out  this  cycle. 

NaN2'''  -f  e'  -*  Na  -1-  N2 
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Figure  26, 


Kinetic  Model  D  —  Og  Minimum 


The  concentrations  of  sodium  neutral  compounds  resulting  from  die  solution  of 
Kinetic  Modd  D. 


46 


The  final  model  is  shown  along  with  all  rate  amstants  in  Table  9.  This  is  the  model  which  will 
be  used  for  the  remainder  of  the  work  and  is  called  Kinetic  Modd  E. 

The  neutral  atomic  sodium  and  the  ions  from  Model  E  are  shown  in  Figure  27.  Sodium  ion 
dominates  above  an  altitude  of  about  IQS  km.  NaN2'*'  b^ins  to  form  as  the  total  atmospheric 
daisity  becomes  substantial  and  mirrors  the  decrease  in  atomic  sodium  with  decreasing  altitude. 
Figure  28  shows  the  sodium  compounds  resulting  from  Model  E.  Excq>t  for  a  decrease  in  the 
top-side  atomic  sodium  doisity,  the  results  are  almost  identical  to  Model  D.  Figure  29  shows 
the  variation  of  sodium  night-glow  with  ozone  doisity.  Here,  the  ozone  doisity  was  varied  by 
multiplication  by  the  factor  0(03. 

The  a=0  limit  has  been  included  in  the  data  in  Figure  29  to  show  the  importance  of  the  reaction 
of  Na  with  O3  in  the  nightglow  process.  Woe  there  no  ozone  at  all,  we  would  see  an  intrasity 
of  around  5  R.  With  low  level  ozrme,  the  intensity  rises  to  about  25  R,  indicating  that  the  ozone 
reaction  is  at  least  a  factor  of  four  more  important  in  producing  the  nightglow  than  is  the  three- 
body  sequence.  Increase  of  the  ozone  by  a  factor  of  ten  increases  the  nigh^ow  by  about  a 
factor  of  five.  The  increase  is  not  linear  becau»  the  process  results  in  an  overall  decrease  of 
fiee  Na.  Also,  increase  in  ozone  increases  the  non-radiative  release  of  NaO  through  reaction 
4.  However,  there  is  still  a  very  strong  correlation  of  the  nightglow  intensity  with  ozone 
density.  In  the  next  section,  we  will  examine  experimental  evid^ice  that  might  support  this 
correlation. 


5.0  SEASONAL  VARIATIONS 


Several  extensive  studies  have  beoi  carried  out  on  the  variations  of  the  sodium  nightglow  at 
various  locations.  One  of  the  most  compressive  was  carried  out  by  Fukuyama  [1977]. 
Among  his  results,  was  the  demonstration  of  a  strong  latitude  dq)endsce  of  the  seasonal 
variations  in  nightglow.  At  mid-latitudes  in  the  northern  hemisphere,  where  we  will  concstrate 
our  own  examination,  he  found  a  maximum  monthly  average  of  about  250  R  in  November  to 
December  and  a  minimum  in  mid-summer  of  about  50  R.  Fukuyama’s  work  is  concerned 
mainly  with  Fourier  analysis  of  the  data  and,  as  such,  does  not  show  the  fine  structure  of  the 
variations  well. 
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Kinetic  Model  E  —  O3  Minimum 


Figure  27.  The  concentation  of  atomic  sodium  and  sodium  ions  resulting  from  die  solution 
of  Kinetic  Model  E. 
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Table  9.  Kinetic  Model  E 


Reaction 


Na  +  O3  -*  NaO  +  O2 


NaO  +  O3  Na02  +  O2 


NaO  +  O3  -*  Na  +  2O2 


NaO  +  HjO  -♦  NaOH  +  OH 


NaO  +  O2  +  N2  Na03  +  N 


NaO  +  CO2  +  N2  NaCOj  +  N 


Na  +  O2  +  N2  Na02  +  N 


Na02  +  O  -*  NaO  +  O2 


Na02  +  H  -*  NaOH  +  O 


NaOH  +  H  Na  +  H2O 


NaOH  +  CO2  +  N2  -*  NaHCOj 


Na02  +  OH  -*  NaOH  +  O2 


Na03  +  O  -*  Na02  +  O2 


NaC03  +  O  -*  Na02  +  CO2 


NaC03  +  H  +  N2  -  NaHC03  +  N2 


NaHC03  +  H  -  Na  +  H2CO3 


Na02  +  CO2  +  N2  “*  NaC03  +  O  +  N 


NaC03  +  Na  +  N2  "*  Na2C03  +  N 


Na2C03  +  O  -♦  NaC03  +  NaO 


Na  +  NO+  -  Na+  +  NO 


7.0(-10) 


1.0(-10) 


2.0(-10) 


2.0(-ll) 


2.0(-10) 


5.3(-30) 


1.3(-27) 


2.0(-30) 


1.0(-13) 


1.9(-14) 


1.9(-28) 


1.0(-16) 


1.4(-9) 


1.0(-9) 


2.5(-31) 


To  tate  a  closer  look  at  the  variations,  we  searched  for  data  from  original  sources. 
Unfortunately,  the  bulk  of  the  nightglow  data  was  recorded  two  or  three  decades  ago  and  much 
of  it  has  been  lost  due  to  storage  problems.  However,  we  were  able  to  locate  about  two  years 
worth  of  data  from  a  few  stations  in  the  early  iwblications  of  World  Data  Center  A  in  Botdd^. 
These  proved  to  be  adequate  for  our  purposes. 

Figure  30  shows  die  combined  data  from  1964  and  196S,  taken  at  Haute  Province,  France. 
Combining  data  is  useful  because  there  are  often  large  gsqis  caused,  presumably,  because  of 
cloudy  conditions  at  the  observation  site.  The  year-to-year  variations  are  actually  quite 
consistent  at  the  same  station.  We  have  also  examined  data  from  Kitt  Peak,  Arizona,  and  find 
that  the  qualitative  agreemoit  is  very  good,  although  latitude  differences  t«id  to  reduce  the 
intensity  of  the  variations  in  the  Kitt  Peak  data  from  those  at  the  45°  north  latitude  of  Haute 
Province. 

The  data  show  two  striking  features  when  processed  with  31-day  running  averages.  First,  there 
are  two  strong  peaks,  one  in  late  March  to  early  April,  and  another  in  early  November.  These 
peaks  are  very  similar  from  year  to  year  and  from  station  to  station.  We  believe  that  these 
peaks,  instead  of  being  caused  by  variations  in  the  ozone,  are  the  results  of  increased  meteoric 
activity.  We  do  not  have  substantiating  evidoice  for  this  belief  at  this  point.  However,  we  note 
that  the  generally  accepted  source  for  the  grains  that  contribute  most  to  the  sodium  in  Earth’s 
atmosphere  is  the  release  of  dust  from  the  evaporation  of  the  ice  in  short  period  comets. 
According  to  Encrenaz,  et  al.  [1989],  Earth  passes  through  the  Lyrids  meteor  swarm  from  comet 
Thatcher  on  10-20  April,  and  the  (Mnids  ftom  Halley,  the  Draconids  from  Giacobini-Zinner, 
and  the  Andromedids  from  Biela  all  in  the  mid-Octoi^r  to  mid-November  time  span.  It  is  of 
interest  that  passage  through  the  Perseids  swarm  from  comet  Swift-Tuttle  is  not  evident  in  Figure 
30.  However,  comet  Swift-TutUe  has  a  relatively  long  period  of  125  years,  and  also  a  perihelion 
somewhat  larger  than  1  A.U.  [Hughes,  1978].  This  may  mean  that  only  the  larger,  and  hence 
brighter  particles  reach  Earth  to  give  rise  to  Ae  Perseids.  We  saw  in  Figure  5  that  the  bulk  of 
the  d^sition  comes  from  particles  with  mass  of  l(-6)  grams.  Further,  the  ratio  of  solar 
radiation  pressure  to  gravitatioiud  attraction  goes  as  the  reciprocal  of  the  particle  radius,  so  that 
smaller  particles  may  be  blown  away,  while  heavier  ones  fall  to  the  sun  to  be  intercepted  by 
Earth.  We  will  leave  off  this  speculation  at  this  point,  concluding  only  that  the  two  pea^  seem 
a  fascinating  subject  for  further  study  from  the  astronomical  viewpoint. 

More  toward  our  purposes  here,  we  note  in  Figure  30  that  there  is  a  strong  sinusoidal  variation 
in  the  nightglow,  neglecting  the  two  prominent  peaks.  The  summer  months  exhibit  a  clear 
decrease  in  intensity  to  about  the  20  R  level  while,  in  winter,  the  intensity  averages  about  140 
R.  This,  we  believe,  is  due  to  changes  in  the  ozone  density  in  the  lower  thermosphere.  In  order 
to  investigate  this  possibility,  we  explored  the  seasonal  variation  in  ozone.  We  began  with 
examining  the  variation  in  total  ozone.  Figure  31  shows  this  variation  as  measured  at  a  station 
near  Haute  Province  at  43°  north.  This  data  has  been  averaged  to  form  monthly  averages  over 
several  years.  As  can  be  seen,  there  is  a  strong  maximum  in  March  and  a  minimum  in  August. 
This  does  not  agree  especially  well  with  the  nightglow  trend  evident  in  Figure  30.  There  are 


52 


Zertth  Intenstiy  (fR) 


Sodium  Night-gbw  at  Haute  Province  (1964-1965) 


Figure  30. 


The  seasonal  variatioD  in  sodium  nightglow  measured  at  Haute  Province,  France, 
for  1964  and  1965. 


53 


Column  Ozone  Density 


amne  good  reasons  for  this.  First,  total  oztxie  column  density  is  measured  in  the  daytime  and 
extrqxdated  to  local  noon.  The  primary  sink  for  02xxie  in  the  mesofi|rfiere  and  thomoqrfiete  in 
the  daytime  is  photo-dissodadtm  via  the  reaction 

O3  +  /n»  -  O2  +  0(^D) , 
while  at  night,  the  reaction  with  atomic  hydrogoi 

O3  +  H  -  OH  +  O2 

dmninates  completdy.  In  the  meso^ete  and  lower  themoq^ere,  the  night-time  increase  in 
Qzoi»  is  som^hing  like  a  foctor  of  10,  so  that  day  measurements  do  not  necessarily  rqnesent 
well  the  night  situaticm  where,  necessarily,  sodium  night-glow  measurements  are  made.  The 
second  reason  why  we  might  not  expect  good  correilation  of  night-glow  witii  total  ozone  is 
simply  the  fact  that  the  total  ozone  is  concentrated  at  lower  altitudes.  The  bulk  of  die  ozone 
colunm  is  in  the  stratosphere  and,  therefore,  the  amount  at  90  km  contributes  negligibly  to  the 
total  column  density. 

To  substantiate  the  lack  of  correlation  of  total  ozone  with  ozone  at  higher  altitudes,  we  cite  the 
work  of  McPeters  [1980]  who  deduced  the  ozone  concoitration  at  1  mbar  pressure  (about  65 
km)  from  Nimbus  4  backscattered  ultraviolet  data.  He  found  a  winter  peak  of  about  8  ftg/g  and 
a  summer  minimum  of  about  4  /tg/g  5(FN,  a  winter  increase  of  a  factor  of  two.  This  is  almost 
exactly  the  opposite  of  the  bdiavior  of  total  ozone,  which  is  a  maximum  in  spring  and  a 
minimum  in 

In  fact,  the  curroit  understanding  of  ozone  density  above  80  km  aiq)ears  at  present  to  be  stetchy 
and  sometimes  contradictory.  The  situation  is  complicated  by  the  existmce  of  a  secondary 
maximum  in  O3  at  the  mesopause,  which  seems  to  be  quite  variable  with  season,  and  even  from 
day  to  day.  There  is  also  a  general  consensus  [see  e.g.  Allen,  et  al. ,  1984],  that  current  models 
do  not  adequately  predict  the  ozone  density  above  80  km  altitude.  Pedi^s  the  most 
comprehoisive  re^ts  to  date  come  from  the  analysis  of  10  years  of  aircraft  02(^A)  emissions 
in  tv^ght,  rqxnted  by  Noxon  [1981].  This  author  summarizes  his  results  as  follows: 

"Betweoi  50  and  70  km  there  appears  to  be  little  variation  (<  ±25%)  whereas 
the  abundance  betweoi  80  and  ^  km  exhibits  a  large  seascmal  change  north  of 
30°N  and  much  less  at  lower  latitude.  At  mid  and  high  latitude  the  colunm 
abundance  above  —80  km  changes  from  <1x10^^  cm*^  in  summer  to  about 
3x10^^  cm’^  in  winter." 

Noxon  made  his  conclusions  by  considering  two  models,  a  "low  ozone"  model  which  fit  the 
decay  in  O2  emissions  (presum^  to  result  from  photolysis  of  O3)  as  a  function  of  z^iith  angle 
well  in  mid-summer,  and  a  "high  ozone"  model  which  fit  well  in  mid-winta.  The  author 
presents  two  curves,  which  we  have  modeled  and  which  are  shown  in  Figure  32.  We  have  used 
these  two  curves  to  model  the  seasonal  variation  as  follows. 
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Altitude  (km) 


Ozone  Density  after  Noxon  [1982] 


Figure  32. 


Miniimim  and  maximum  values  of  ozone  density  in  the  mesosphere  and  Iowa* 
troposphere,  modded  afta  the  results  of  Noxon  [1982]. 
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where  m  is  the  month.  This  model,  when  used  to  rqilace  the  O3  density  in  Model  E,  gives  rise 
to  the  nightglow  variaticHi  shown  in  Figure  33.  The  variation  is  in  remarkably  good  agreement 
with  the  variation  in  nightglow  at  Haute  Province,  shown  in  Figure  30,  save  again  for  the  two 
peaks  we  believe  to  be  due  to  m^eor  showers.  The  winter  maximum  is  about  130  R,  compared 
to  perlu^  140  R  in  the  Haute  Province  data,  and  the  summer  minimum  is  around  30  R  in  the 
model,  compared  to  around  40  R  in  the  data. 

In  Figure  34,  we  show  the  emission  profile  as  a  function  of  season,  and  Figure  35  shows  the 
variation  in  Ae  sodium  layer  caused  by  the  variation  in  We  see  that  the  nuxiel  predicts  a 
slight  lowering  of  the  layer  in  the  summer  with  a  concurrent  increase  in  sodium  density.  The 
increase  is  actually  quite  slight,  amounting  to  about  20%  in  terms  of  column  density.  Megie  and 
Blamount  rqxnt  that  the  sodium  layer  at  Haute  Province  is  actually  quite  ctmstant,  excqH  for 
a  rather  sharp  increase  in  late  November.  This  increase  would  not  seem  to  correlate  wiA  our 
model  results,  however,  it  could  be  due  to  the  increase  in  meteoric  activity  since  it  does  appear 
to  correlate  with  the  second  anomalous  peak  in  the  nightglow  data.  The  observed  increase  in 
free  Na  is  more  than  a  factor  of  10,  which  is  far  too  much  to  be  explained  by  changes  in  the 
[O3]. 

There  are  at  least  two  variables  which  we  have  omitted  in  this  examination.  One  of  these  is 
temperature.  The  variation  in  temperature  will  affect  the  rates  of  all  kinetic  reactions  and  will 
also  influence  the  thermal  diffusion  term  in  the  diffusion  equation.  We  do  not  rule  out  the 
possible  influence  of  temperature  effects,  however,  we  believe  from  simple  calculations  that  they 
would  be  of  secondary  importance.  Finally,  the  nightglow  and  the  sodium  density  are  both 
strongly  influoiced  by  the  atomic  oxygoi  concentraticm.  This  is  simply  a  result  of  the  fiict  that 
the  reaction  between  Na  and  O3  is  the  rate  limiting  stq>  in  the  Chapman  process.  The  atomic 
hydrogen  doisity  is  also  of  potoitial  importance  in  rqploiishing  frre  sodium.  Unfortunately, 
neither  of  these  is  a  particularly  well  known  quantity.  We  believe  that  in  view  of  the  uncertainty 
in  the  rate  coefficients  used,  the  assumption  of  no  seasonal  variation  in  the  atomic  atmospheric 
species  is  reasonable.  Also,  the  COSPAR  International  Reference  Atmosphere  of  1986,  rqmrts 
no  seasonal  variation  of  atomic  oxygen  at  90  km. 
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Figure  35.  The  seasonal  variation  in  free  sodium  using  the  modeled  Qj  profiles. 
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6.0  SUMMARY 


We  have  developed  a  comprehoisive  model  of  atmosj^eiic  sodium  deposition,  b^inning  only 
with  assumptions  of  the  total  mass  dqx>sition  rate.  The  model  yields  a  sodium  layer  that  is  in 
good  agreemoit  with  experimoital  results.  We  have  used  the  mo^  to  predict  sodium  nightglow 
intensity  and  have  succeeded  in  duplicating  the  measured  seasonal  variation  in  that  intensity  by 
varying  the  lower  tropospheric  ozone  in  accordance  with  experimentally  determined  ozone 
minima  and  maxima.  The  results  agree  quite  well  with  the  measured  ni^tglow  at  the  4S‘’N 
station  we  examined. 

The  model  itsdf  is  pleasing  in  that  it  is  based  to  a  large  extent  on  first  principles,  b^inning 
from  an  assumed  and  relatively  well  known  mass  input  rate  of  cosmic  dust  to  Earth’s 
atmosphere.  This  modd  is  equally  valid  for  any  other  memllic  spedes  that  arises  from  meteoric 
particle  ablation  and  is  valid  without  modification  at  any  altitude.  We  intend  to  use  the  basic 
model  to  study  the  deposition  of  other  metals  and  espedally  of  ions  at  higher  altitudes.  As  is 
generally  the  case,  the  introduction  of  a  comprdi«isive  kinetic  scheme  has  required  some 
assumptions  about  various  rate  constants.  However,  in  this  present  model,  we  have  obtained 
a  scheme  which  is  not  strongly  dq)endent  on  any  of  the  predse  rates  for  the  “closure” 
reactions,  the  rates  for  which  are  not  well  known.  Tlie  basic  results  of  this  model  would  remain 
the  same  with  variation  of  rate  constants  within  reasonable  bounds. 

Finally,  the  good  correlation  of  modded  sodium  nightglow  with  experimental  results  strongly 
suppo^  the  conclusion  that  lower  thermospheric  ozone  abundance  is  the  driving  force  in  the 
seasonal  variation  in  sodium  nightglow.  This  result  is  important  in  that  this  region  of  the 
atmosphere  is  not  well  understood  goierally,  and  current  modds  do  not  adequately  predict  the 
chemistry  there.  We  believe  that  models  such  as  this  one,  which  attempt  to  identify  causative 
agents  for  well  known  behavior,  can  contribute  to  a  comprehoisive  understanding  of  the  region. 
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APPENDK.  KINETIC  MODEL  SOLUTION  METHOD 


The  problem  of  solving  the  steady  state  kinetic  modds  rapidly  becomes  complicated  with  an 
increasing  number  of  reactions  and  constituents.  We  have  written  a  program  called  KINO  which 
assists  in  this  task.  The  routine  takes  as  input  a  file  describing  the  kinetic  model.  The  file  for 
the  most  complicated  model  used  in  this  investigation,  Model  E,  is  shown  below. 


SODIUM  MODEL  E 

01) 

Na  +  03  ^  NaO  +  02 

7.0E-10 

02) 

NaO  +  0  Na  +  02 

$ 

1.0E-10 

03) 

NaO  +  03  =  Na02  +  02 

$ 

2.0E-10 

04) 

NaO  +  03  =  Na  +  02  -i-  02 

$ 

2.0E-11 

05) 

NaO  +  H20  =  NaOH  +  OH 

$ 

2.0E-10 

06) 

NaO  +  02  +  N2  =  Na03  +  N2 

$ 

5.3E-30 

07) 

NaO  +  C02  +  N2  =  NaC03  +  N2 

$ 

1.3E-27 

08) 

Na  ^  02  +  N2  ^  Na02  +  N2 

$ 

2.0E-30 

09) 

Na  +  02+  =  A/a+  +  02 

$ 

1.4E-9 

10) 

Na  +  /V0+  =  /Va+  +  NO 

$ 

1.0e-9 

11) 

Na+  +  N2  +  N2  =  NaN2+  +  N2 

4 

2.5e-31 

12) 

NaN2+  +  e-  -  Na  +  N2 

$ 

1.0e-6 

13) 

Na02  +  0  =  NaO  +  02 

$ 

1.0E-13 

14) 

Na02  +  =  NaOH  +  0 

S 

1.9E-14 

15) 

NaOH  +  H  =  Na  +  H20 

S 

1.0e-11 

16) 

NaOH  +  002  +  /V2  =  NaHC03  +  N2 

$ 

1.9E-28 

17) 

Na02  +  OH  =  NaOH  +  02 

$ 

1.0E-10 

18) 

Na03  +  0  =  Na  +  02  +  02 

$ 

2.0E-14 

19) 

Na03  +  0  =  Na02  +  02 

$ 

2.0E-14 

20) 

NaC03  +  0  =  Na02  +  002 

$ 

1.0E-11 

21) 

NaC03  +  H  +  N2  =  NaHC03  +  N2 

$ 

1.0E-28 

22) 

NaHC03  +  //  =  /Va  +  H2C03 

$ 

1.0E-11 

23) 

Na02  +  002  +  N2  =  NaC03  +  0  +  N2 

$ 

1.0E-30 

24) 

NaC03  +  Na  +  N2  =  Na2C03  +  N2 

1.0E-28 

25) 

Na2C03  +  0  =  NaO  +  NaC03 

$ 

1.0E-16 

The  input  file  contains  a  line  consisting  of  the  reactants  and  products,  which  are  sq)arated  by 
an  '=’  sign.  After  these,  a  *$’  sign  indicates  that  the  next  value  is  the  rate  coefficient.  The 
KINO  program  reads  this  file  and  tallies  up  all  reactants  and  products  along  with  the  appropriate 
rate  constant.  The  program  is  peihaps  b^t  described  by  showing  its  various  output. 

The  first  ou^ut  is  a  list  of  all  different  species  encountered  in  the  input  file.  Along  with  this, 
is  an  indication  of  which  reactions  these  species  are  involved  in.  The  list  below  corresponds  to 
the  input  file  above.  A  second  table  (not  shown  in  this  output)  keeps  tracks  of  how  many  of 
each  species  are  created  and/or  destroyed  by  each  reaction. 
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After  presenting  this  table,  the  program  can  do  a  mass  check  and  charge  check  to  make  sure 
there  have  beoi  no  typos  in  the  input  file.  The  results  are  given  below.  Charge  is  indicated  in 
the  file  by  a  ‘+’  or  a  following  the  species  in  the  input  file. 


FINISHED  WITH  INPUT 

DO  YOU  WANT  A  CHARGE  BALANCE  CHECK?!  1  =  YES! 


CHECKING  CHARGE  BALANCE . 


SPECIES 

1 

Na 

CHARGE 

0 

SPECIES 

2 

03 

CHARGE 

0 

SPECIES 

3 

NaO 

CHARGE 

0 

SPECIES 

4 

02 

CHARGE 

0 

SPECIES 

5 

0 

CHARGE 

0 

SPECIES 

6 

Na02 

CHARGE 

0 

SPECIES 

7 

H20 

CHARGE 

0 

SPECIES 

8 

NaOH 

CHARGE 

0 

SPECIES 

9 

OH 

CHARGE 

0 

SPECIES 

10 

N2 

CHARGE 

0 

SPECIES 

11 

Na03 

CHARGE 

0 

SPECIES 

12 

C02 

CHARGE 

0 

SPECIES 

13 

NaC03 

CHARGE 

0 

SPECIES 

14 

02+ 

CHARGE 

1 

SPECIES 

15 

Na  + 

CHARGE 

1 

SPECIES 

16 

N0+ 

CHARGE 

1 

SPECIES 

17 

NO 

CHARGE 

0 

SPECIES 

18 

NaN2+ 

CHARGE 

1 

SPECIES 

19 

e- 

CHARGE 

-1 

SPECIES 

20 

H 

CHARGE 

0 

SPECIES 

21 

NaHC03 

CHARGE 

0 
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spsaes  22  H2C03  CHARGE  0 

SPECmS  23  Na2C03  CHARGE  0 

RXN.  1  /V«  +  03  «  NaO  +  02 

RXN.  2  NaO  +  O  »  Ato  +  02 

RXN.  3  NaO  -t-  03  -  Na02  -t-  02 

RXN.  4  NaO  +  03  »  +  02  +  02 

RXN.  S  NaO  +  H20  «  NaOH  +  OH 

RXN.  6  NaO  02  N2  ^  Na03  N2 

RXN.  7  NaO  +  002  ■¥  N2  ^  NaC03  +  N2 

RXN.  a  Ato  +  02  +  iV2  -  Na02  +  N2 

RXN.  9  Na  +  02+  »  Na+  +  02 

RXN.  10  Na  +  N0+  »  Na+  +  NO 

RXN.  11  Na+  +  N2  +  N2  ^  NaN2+  -i-  N2 

RXN.  12  NaN2+  +  e-  •  Na  +  N2 

RXN.  13  Na02  +  O  »  NaO  +  02 

RXN.  14  Na02  +  H  =  NaOH  +  O 

RXN.  IS  NaOH  +  H  =  Na  +  H20 

RXN.  16  NaOH  +  002  +  N2  ^  NaHC03  +  N2 

RXN.  17  Na02  +  OH  ^  NaOH  +  02 

RXN.  18  Na03  +  O  ^  Na  +  02  +  02 

RXN.  19  Na03  +  O  =  Na02  +  02 

RXN.  20  NaC03  +  O  «  Na02  +  C02 

RXN.  21  NaC03  +  H  +  N2  =  NaHC03  +  N2 

RXN.  22  NaHC03  +  H  =  Na  +  H2C03 

RXN.  23  Na02  +  C02  +  /V2  =  NaC03  +  O  +  N2 

RXN.  24  NaC03  +  Na  +  N2  ^  Na2C03  +  N2 

RXN.  25  Na2C03  +  O  =  NaO  +  NaC03 


CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 
CHARGE  0 


DO  YOU  WANT  A  MASS  BALANC  ^CK?(1  =  YES1 


Sp. 

1 

Na 

Mass 

22 

Sp. 

2 

03 

Mass 

48 

Sp. 

3 

NaO 

Mass 

38 

Sp. 

4 

02 

Mass 

32 

Sp. 

5 

0 

Mass 

16 

Sp. 

6 

Na02 

Mass 

54 

Sp. 

7 

H20 

Mass 

18 

Sp. 

8 

NaOH 

Mass 

39 

9 

OH 

Mass 

17 

Sp. 

10 

N2 

Mass 

28 

Sp. 

11 

Na03 

Mass 

70 

Sp. 

12 

C02 

Mass 

44 

Sp. 

13 

NaC03 

Mass 

82 

14 

02+ 

Mass 

32 
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Sp.  IS 

Na+ 

Mass 

22 

Sp.  16 

N0+ 

Mass 

30 

Sp.  17 

NO 

Mass 

30 

Sp.  18 

NaN2+ 

Mass 

50 

Sp.  19 

a- 

Mass 

0 

Sp.  20 

H 

Mass 

1 

Sp.  21 

NaHC03 

Mass 

83 

Sp.  22 

H2C03 

Mass 

62 

Sp.  23 

Na2C03 

Mass 

104 

RXN.  1 

Na  +  03 

=  NaO  +  02 

RXN.  2 

NaO  +  0 

^  Na  +  02 

RXN.  3  NaO  +  03  «  Na02  -i-  02 

RXN.  4  NaO  +  03  ^  Na  +  02  +  02 

RXN.  5  NaO  +  H20  »  NaOH  +  OH 

RXN.  6  NaO  +  02  +  N2  ^  Na03  +  N2 

RXN.  7  NaO  +  C02  +  N2  NaC03  +  N2 

RXN.  8  Na  +  02  +  N2  ^  Na02  +  N2 

RXN.  9  Na  +  02+  =  Na+  +  02 

RXN.  10  Na  +  NO+  =  Na+  +  NO 

RXN.  11  Na+  +  N2  +  N2  =  NaN2+  +  N2 

RXN.  12  NaN2+  +  e-  ^  Na  +  N2 

RXN.  13  Na02  +  O  «  NaO  +  02 

RXN.  14  Na02  +  //  »  NaOH  +  O 

RXN.  IS  NaOH  +  H  ^  Na  +  H20 

RXN.  16  NaOH  +  C02  +  N2  ^  NaHC03  +  N2 

RXN.  17  Na02  +  OH  ^  NaOH  +  02 

RXN.  18  Na03  +  O  =  Na  +  02  +  02 

RXN.  19  Na03  +  O  «  Na02  +  02 

RXN.  20  NaC03  +  O  =  Na02  +  C02 

RXN.  21  NaC03  +  H  +  N2  =  NaHC03  +  N2 

RXN.  22  NaHC03  +  H  =  Na  +  H2C03 

RXN.  23  Na02  +  C02  +  N2  ^  NaC03  +  O  +  N2 

RXN.  24  NaC03  +  Na  +  N2  ^  Na2C03  +  N2 

RXN.  25  Na2C03  +  O  =  NaO  +  NaC03 


MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 

MASS  0 


The  mass  is  determined  by  adding  up  all  species  in  each  reaction,  and  a  zero  total  mass  for  each 
reaction  indicates  a  valid  reaction,  as  does  a  zero  total  charge.  After  this  output,  the  routine  can 
write  either  a  time-dependrat  subroutine  using  Runge-Kutta  or  a  steady  state  routine.  For  the 
steady  state  routine,  a  matrix  is  goierated.  To  understand  the  matrix  cmicept,  consider  the 
following  simple  reaction  scheme. 
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Na  *  NaO  *  02 
NaO  *  0-^Na  *  02 

In  tiiis  sinqde  sdheme,  Acre  are  two  NaX  q)ecies,  Na  and  NaO.  The  steady  state  equations  can 
dius  be  written 


-  0  -  -ikiiAfallOjl  +  ki[NaO][0] 
-  0  -  tjIMilIOal  -  k^lNaOnOl 


AldKNigh  diis  case  is  trivial  to  solve,  we  see  that  the  problem  can  be  written  in  die  form  of 
linear  system  in  whidi  die  elements  of  an  array  a^  are  the  sums  of  all  reaction  rates  that  remove 
die  conqwnent  i  from  the  system,  and  the  a^  are  the  sums  of  all  those  that  change  component 
j  into  component  i.  We  can  also  see  that  die  system  is  overdetermined,  since  we  can  multiply 
all  the  components  by  a  scaler  and  sdll  have  a  solution. 


\\  <*12  <*13* 

■«r 

<^1  <*22  <*23 

**2 

_<*31  <*32  <*33_ 

/*3. 

[01 


The  equation  above  rqnesents  a  three-component  steady  state  matrix.  To  solve  the  problem, 
we  tate  the  value  of  one  of  the  compcments,  nj,  to  be  unity.  Then,  we  solve  the  n-1  linear 
equation  to  find  die  components  02,  n3,....nn.  Finally,  knowing  the  totil  of  all  components,  we 
can  normalize  the  soluticm.  The  KINO  program  writes  a  subroutine  whidi  fills  up  a  matrix, 
such  as  the  one  above.  A  sample  subroutine  correspcmding  to  Model  E  is  givoi  bdow.  Li  that 
routine,  die  qiecies  that  must  be  qiecified  externally,  Le.  the  constant  components  of  die 
atmaqdiere  in  this  case,  are  listed  and  must  be  supplied  to  the  routine  through  commmi 
statements.  The  KINO  program  asks  the  user  for  a  string,  in  diis  case  ’Na’,  by  whidi  it  may 
differentiate  omstant  axistituoits  fiom  variable  constituents. 


SUBROUTINE  FILLMA  TfAMA  V 
C 

IMPUCrr  REAL*8  (A-H,0-Z) 

C 

DIMENSION  AMA  T(10, 10) 

C 
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« 


C  THE  FOUOymQ  ARE  VARIABLES  IN  THE  PROBLEM 
C 

C  1  2Nm 
C  2  3NmO 
C  3  4  Nm02 

C  4  SNmOH 
C  SB  Ns03 
C  6  7  NaC03 

C  7  8  Nap 

C  8  B  NaN2p 
C  9  10  NaHC03 

C  10  11  Na2C03 

C 

C  THE  FOLLOWING  DENSITIES  MUST  BE  SUPPLIED 
C 

C  03 
C  02 
C  O 
C  H20 
C  OH 
C  N2 
C  C02 
C  02p 
C  NOp 
C  am 
C  H 
C 

C  THE  FOLLOWING  REACTIONS  ARE  INCLUDED 
C 

C  1.  Na  ^  03  ^  NaO  -i-  02 
C  2.  NaO  -I-  O  =  yVa  +  02 
C  3.  NaO  +  03  =  Na02  +  02 
C  4.  NaO  +  03  ^  Na  +  02  +  02 
C  5.  NaO  +  H20  «  NaOH  +  OH 
C  6.  NaO  +  02  +  N2  =  Na03  +  N2 
C  7.  NaO  +  C02  +  N2  =  NaC03  +  N2 
C  8.  Na  +  02  +  N2  =  Na02  +  N2 
C  9.  Na  +  02+  =  Na+  +  02 
C  10.  Na  +  N0+  =  Na+  +  NO 
C  11.  Na+  +  N2  +  N2  =  NaN2+  +  N2 
C  12.  NaN2+  +  e-  ^  Na  +  N2 
C  13.  Na02  +  O  =  NaO  +  02 
C  14.  Na02  +  //  =  NaOH  +  O 
C  15.  NaOH  +  H  ^  Na  +  H20 
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C  16.  NsOH  ■¥  C02  +  N2  m  NaHCOB  +  N2 
C  17.  Nm02  +  oh  ~  N»OH  +  02 
C  18.  Ns03  O  Hs  +  02  +  02 
C  19.  Ha03  +  O  •  N»02  +  02 
C  20.  NaC03  +  O  -  Af«02  +  C02 
C  21.  NaC03  ^  H  ^  N2  ^  NmHC03  +  N2 
C  22.  NmHC03  +  //«/»•  +  H2C03 
C  23.  Na02  +  C02  +  Af2  -  NaC03  +  O  +  iV2 
C  24.  NaC03  ^  Na  N2  ^  Na2C03  •(-  N2 
C  25.  Na2C03  +  O  »  NaO  ^  NaC03 
C 

REALMS  03 
REALMS  02 
REALMS  O 
REAL*8H20 
REAL*8  0H 
REAL*8N2 
REAL*8C02 
REAL*8  02p 
REAL*8NOp 
REAL*8am 
REAL*8H 
C 

03-=  003 
02-002 
O-DO 
H20-DH20 
OH-DOH 
N2-DN2 
C02-DC02 
02p=D02p 
NOp-DNOp 
am— Dam 
H-DH 
C 

R01-  7.0000E-10 
R02-  1.0000E-10 
R03-  2.0000E-10 
R04-  2.0000E-11 
R05=  2.0000E-10 
R06-  5.3000E-30 
R07-  1.3000E-27 
R08-  2.0000E-30 
R09-  1.4000E-09 


70 


RIO*  1.0000£-09 
R11*  2.S000E-31 
R12*  1.0000E-06 
R13*  1.0000E-13 
R14*  1.9000£-14 
R1S*  1.0000E-11 
R16*  1.9000E-28 
R17*  1.0000E-10 
R18*  2.0000E-14 
R19*  2.0000E-14 
R20*  1.0000E-11 
R21*  1.0000E-28 
R22*  1.0000E-11 
R23*  1.0000E-30 
R24*  1.0000E-28 
R25*  1.0000E-16 
C 

AMATf01,01i*-R01*03‘R08*02^N2-R09*02p-R10*NOp-R24*N2 

AMATf01.02)*  ■¥R02*0-¥R04*03 

AMAT(0h03}*0.0 

AMA  Tf0 1,04)*+  R1S*H 

AMAT(01,0S)*  +R18*0 

AMAT(01,06)*0.0 

AMATf01,07)  =  0.0 

AMAT(01,08)*  +R12*0m 

AMA  T(01,09)  =  +  R22*H 

AMAT(01.10)*0.0 

AMAT(02,01)*  +R01*03 

AMArf02,02)*-R02*0-R03*03-R04*03-ROS*H20-R06*02*N2-R07*N2*C02 

AMA  T(02,03)  *  +  R13*0 

AMAT(02,04)*0.0 

AMAT(02,0S)*0.0 

AMAT(02,06)*0.0 

AMATf02,07)*0.0 

AMAT(02,08)*0.0 

AMAT(02,09)  =  0.0 

AMA  7(02, 10)  *  +  R25*0 

AMA  7(03,01)  =  +  R08*02*N2 

AMA 7(03,02)  *  +  R03*03 

AMA  7(03,03)  =  -R13*0-R14*H-R17*OH-R23*N2*C02 

AMA7(03,04)*0.0 

AMA 7(03,05)  *+R19*0 

AMA7(03,08)*  +R20*0 

AMA7(03,07)*0.0 
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AMATi03,08)^0.0 

AMATf03,09)»0.0 

AMATfa3,Wm0.0 

AMATf04,01)^0.0 

AMATm,02)^  ^R0S*H20 

AMAT(04,03)^  ■¥R14*H-^R17*OH 

AMA 7(04,04)  «  -R1S*H-R16*N2*C02 

AMAT(04,0S)^0.0 

AMAT(04,06)^0.0 

AMATf04,07)^0.0 

AMAT(04,08}^0.0 

AMAT(04,09}^0.0 

AlliAT(04,10)»0.0 

AMAT(0S,01)^0.0 

AMAT(0S,02)^  ■¥R06*02^N2 

AMAT(0S,03)^0.0 

AMAT(0S,04J  =  0.0 

AMA  TfOS,OS)  >  -R18*0-R19*0 

AMATf0S,06J  =  0.0 

AMAT(0S,07)  =  0.0 

AMAT(0S,08)^0.0 

AMAT(0S,09)^0.0 

AMAT(0S,10)^0.0 

AMAT(06,01)  =  0.0 

AMA 7(06,02}  *  +  R07*N2*C02 

AMA  7(06,03)  =  +  R23*N2*C02 

AMA7(06,04)  =  0.0 

AMA7(06,05)  =  0,0 

AMA 7(06,06)  »  -R20*0-R21*N2*H-R24*N2 

AMA7(06,07)  =  0.0 

AMA7(06,08)  =  0.0 

AMA7(06,09)  =  0.0 

AMA  7(06, 10)  »  -i-  R25*0 

AMA7(07,01)=  •¥R09*02p-^R10*NOp 

AMA7(07,02)  =0.0 

AMA7(07,03)  =  0.0 

AMA7(07,04)  =  0.0 

AMA7(07,05J  =  0.0 

AMA7(07,06)^0.0 

AMA  7(07,07)  =  -R1 1  •N2*N2 

AMA7(07,08)  =0.0 

AMA7(07,09)  =  0.0 

AMA7(07,10J^0.0 

AMA7(08,01)  =  0.0 
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AMAT(08,02)^0.0 
AMATfOa^Oaj^rO.O 
AMAT(08,<W^0.0 
AMArf08,0S)^0.0 
AMATf08,06)^0.0 
AMA  T(0B.07)  »  +  /?7  7  *N2*N2 
AMA  7(08,08}  -  ‘M2*mn 
AMAT(08,09}^0.0 
AMAT(08,10}^0.0 
AMAT(09,01}^0.0 
AlillAT(09,02}^0.0 
AMAT(09,03}^0.0 
AMA 7(09,04}  =  +  R16*N2*C02 
AMA7(09,0S}  =  0.0 
AMA 7(09,06}  ^•*■821  *N2*H 
AMA7(09,07}^0.0 
AMA7(09,08}^0.0 
AMA 7(09,09}  »  -R22*H 
AMA7(09,10}^0.0 
AMA  7(10,01}  =  +  R24*N2 
AMA7(  10,02}  =  0.0 
AMA7(10,03}  =  0.0 
AMA7(10,04}  =  0.0 
AMA7(10,05}  =  0.0 
AMA7(10,0€}=  •k-R24*N2 
AMA7(  10,07}  =  0.0 
AMA7(10,08}  =  0.0 
AMA7(10,09}  =  0.0 
AMA  7(10, 10}  =  -R25*0 
C 

RETURN 

END 
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